{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _*Credit Risk Analysis*_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "This tutorial shows how quantum algorithms can be used for credit risk analysis.\n",
    "More precisely, how Quantum Amplitude Estimation (QAE) can be used to estimate risk measures with a quadratic speed-up over classical Monte Carlo simulation.\n",
    "The tutorial is based on the following papers:\n",
    "\n",
    "- [Quantum Risk Analysis. Stefan Woerner, Daniel J. Egger.](https://www.nature.com/articles/s41534-019-0130-6) [Woerner2019]\n",
    "- [Credit Risk Analysis using Quantum Computers. Egger et al. (2019)](https://arxiv.org/abs/1907.03044) [Egger2019]\n",
    "\n",
    "A general introduction to QAE can be found in the following paper:\n",
    "\n",
    "- [Quantum Amplitude Amplification and Estimation. Gilles Brassard et al.](http://arxiv.org/abs/quant-ph/0005055)\n",
    "\n",
    "The structure of the tutorial is as follows:\n",
    "\n",
    "1. [Problem Definition](#Problem-Definition)\n",
    "2. [Uncertainty Model](#Uncertainty-Model)\n",
    "3. [Expected Loss](#Expected-Loss)\n",
    "4. [Cumulative Distribution Function](#Cumulative-Distribution-Function)\n",
    "5. [Value at Risk](#Value-at-Risk)\n",
    "6. [Conditional Value at Risk](#Conditional-Value-at-Risk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from qiskit import QuantumRegister, QuantumCircuit, Aer, execute\n",
    "from qiskit.circuit.library import IntegerComparator\n",
    "from qiskit.aqua.algorithms import IterativeAmplitudeEstimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Definition\n",
    "\n",
    "In this tutorial we want to analyze the credit risk of a portfolio of $K$ assets.\n",
    "The default probability of every asset $k$ follows a *Gaussian Conditional Independence* model, i.e., given a value $z$ sampled from a latent random variable $Z$ following a standard normal distribution, the default probability of asset $k$ is given by\n",
    "\n",
    "$$p_k(z) = F\\left( \\frac{F^{-1}(p_k^0) - \\sqrt{\\rho_k}z}{\\sqrt{1 - \\rho_k}} \\right) $$\n",
    "\n",
    "where $F$ denotes the cumulative distribution function of $Z$, $p_k^0$ is the default probability of asset $k$ for $z=0$ and $\\rho_k$ is the sensitivity of the default probability of asset $k$ with respect to $Z$. Thus, given a concrete realization of $Z$ the individual default events are assumed to be independent from each other.\n",
    "\n",
    "We are interested in analyzing risk measures of the total loss\n",
    "\n",
    "$$ L = \\sum_{k=1}^K \\lambda_k X_k(Z) $$\n",
    "\n",
    "where $\\lambda_k$ denotes the _loss given default_ of asset $k$, and given $Z$, $X_k(Z)$ denotes a Bernoulli variable representing the default event of asset $k$. More precisely, we are interested in the expected value $\\mathbb{E}[L]$, the Value at Risk (VaR) of $L$ and the Conditional Value at Risk of $L$ (also called Expected Shortfall). Where VaR and CVaR are defined as\n",
    "\n",
    "$$ \\text{VaR}_{\\alpha}(L) = \\inf \\{ x \\mid \\mathbb{P}[L <= x] \\geq 1 - \\alpha \\}$$\n",
    "\n",
    "with confidence level $\\alpha \\in [0, 1]$, and\n",
    "\n",
    "$$ \\text{CVaR}_{\\alpha}(L) = \\mathbb{E}[ L \\mid L \\geq \\text{VaR}_{\\alpha}(L) ].$$\n",
    "\n",
    "For more details on the considered model, see, e.g.,<br>\n",
    "<a href=\"https://arxiv.org/abs/1412.1183\">Regulatory Capital Modeling for Credit Risk. Marek Rutkowski, Silvio Tarca</a>\n",
    "\n",
    "\n",
    "\n",
    "The problem is defined by the following parameters:\n",
    "- number of qubits used to represent $Z$, denoted by $n_z$\n",
    "- truncation value for $Z$, denoted by $z_{\\text{max}}$, i.e., Z is assumed to take $2^{n_z}$ equidistant values in $\\{-z_{max}, ..., +z_{max}\\}$ \n",
    "- the base default probabilities for each asset $p_0^k \\in (0, 1)$, $k=1, ..., K$\n",
    "- sensitivities of the default probabilities with respect to $Z$, denoted by $\\rho_k \\in [0, 1)$\n",
    "- loss given default for asset $k$, denoted by $\\lambda_k$\n",
    "- confidence level for VaR / CVaR $\\alpha \\in [0, 1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set problem parameters\n",
    "n_z = 2\n",
    "z_max = 2\n",
    "z_values = np.linspace(-z_max, z_max, 2**n_z)\n",
    "p_zeros = [0.15, 0.25]\n",
    "rhos = [0.1, 0.05]\n",
    "lgd = [1, 2]\n",
    "K = len(p_zeros)\n",
    "alpha = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We now construct a circuit that loads the uncertainty model. This can be achieved by creating a quantum state in a register of $n_z$ qubits that represents $Z$ following a standard normal distribution. This state is then used to control single qubit Y-rotations on a second qubit register of $K$ qubits, where a $|1\\rangle$ state of qubit $k$ represents the default event of asset $k$. The resulting quantum state can be written as\n",
    "\n",
    "$$ |\\Psi\\rangle = \\sum_{i=0}^{2^{n_z}-1} \\sqrt{p_z^i} |z_i \\rangle \\bigotimes_{k=1}^K \n",
    "\\left( \\sqrt{1 - p_k(z_i)}|0\\rangle + \\sqrt{p_k(z_i)}|1\\rangle\\right),$$\n",
    "\n",
    "where we denote by $z_i$ the $i$-th value of the discretized and truncated $Z$ [Egger2019]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.finance.applications import GaussianConditionalIndependenceModel as GCI\n",
    "u = GCI(n_z, z_max, p_zeros, rhos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌─────────┐┌─────────┐\n",
       "q_0: ┤0      ├┤0        ├┤0        ├\n",
       "     │  P(X) ││         ││         │\n",
       "q_1: ┤1      ├┤1 LinRot ├┤1        ├\n",
       "     └───────┘│         ││  LinRot │\n",
       "q_2: ─────────┤2        ├┤         ├\n",
       "              └─────────┘│         │\n",
       "q_3: ────────────────────┤2        ├\n",
       "                         └─────────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌─────────┐┌─────────┐\n",
       "q_0: ┤0      ├┤0        ├┤0        ├\n",
       "     │  P(X) ││         ││         │\n",
       "q_1: ┤1      ├┤1 LinRot ├┤1        ├\n",
       "     └───────┘│         ││  LinRot │\n",
       "q_2: ─────────┤2        ├┤         ├\n",
       "              └─────────┘│         │\n",
       "q_3: ────────────────────┤2        ├\n",
       "                         └─────────┘"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the simulator to validate the circuit that constructs $|\\Psi\\rangle$ and compute the corresponding exact values for\n",
    "- expected loss $\\mathbb{E}[L]$\n",
    "- PDF and CDF of $L$ \n",
    "- value at risk $VaR(L)$ and corresponding probability\n",
    "- conditional value at risk $CVaR(L)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the circuit and analyze the results\n",
    "job = execute(u, backend=Aer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analyze uncertainty circuit and determine exact solutions\n",
    "p_z = np.zeros(2**n_z)\n",
    "p_default = np.zeros(K)\n",
    "values = []\n",
    "probabilities = []\n",
    "num_qubits = u.num_qubits\n",
    "for i, a in enumerate(job.result().get_statevector()):\n",
    "    \n",
    "    # get binary representation\n",
    "    b = ('{0:0%sb}' % num_qubits).format(i)\n",
    "    prob = np.abs(a)**2\n",
    "\n",
    "    # extract value of Z and corresponding probability    \n",
    "    i_normal = int(b[-n_z:], 2)\n",
    "    p_z[i_normal] += prob\n",
    "\n",
    "    # determine overall default probability for k \n",
    "    loss = 0\n",
    "    for k in range(K):\n",
    "        if b[K - k - 1] == '1':\n",
    "            p_default[k] += prob\n",
    "            loss += lgd[k]\n",
    "    values += [loss]\n",
    "    probabilities += [prob]   \n",
    "\n",
    "values = np.array(values)\n",
    "probabilities = np.array(probabilities)\n",
    "    \n",
    "expected_loss = np.dot(values, probabilities)\n",
    "\n",
    "losses = np.sort(np.unique(values))\n",
    "pdf = np.zeros(len(losses))\n",
    "for i, v in enumerate(losses):\n",
    "    pdf[i] += sum(probabilities[values == v])\n",
    "cdf = np.cumsum(pdf)\n",
    "\n",
    "i_var = np.argmax(cdf >= 1-alpha)\n",
    "exact_var = losses[i_var]\n",
    "exact_cvar = np.dot(pdf[(i_var+1):], losses[(i_var+1):])/sum(pdf[(i_var+1):])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected Loss E[L]:                0.6409\n",
      "Value at Risk VaR[L]:              2.0000\n",
      "P[L <= VaR[L]]:                    0.9591\n",
      "Conditional Value at Risk CVaR[L]: 3.0000\n"
     ]
    }
   ],
   "source": [
    "print('Expected Loss E[L]:                %.4f' % expected_loss)\n",
    "print('Value at Risk VaR[L]:              %.4f' % exact_var)\n",
    "print('P[L <= VaR[L]]:                    %.4f' % cdf[exact_var])\n",
    "print('Conditional Value at Risk CVaR[L]: %.4f' % exact_cvar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAElCAYAAAAhjw8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZgUxfnA8e/LKct9iAqCiAqoqESRAB6ggBGFEBGBVVG8iBo03tcPFY8YD0CTeEJMQDwWD4hREBSS9UINoCRBDjWKHIIIcu+CHO/vj+qBYXaOnrN3d97P8/QzO9VV3e92z05td3VViapijDHGpKNK0AEYY4yp+KwyMcYYkzarTIwxxqTNKhNjjDFps8rEGGNM2qwyMcYYkzarTIxJg4iMFxEVkVYB7X+kt//uEekqIsVBxBQWQ6DHxuSWVSYmZd4XRYXvqCQiS0O/i7fsEJF1IvJfEZkoIueJSI0s7bu4oh7DWBWZyU/Vgg7AmHLkD8AG3D9Z9YC2wDnAhcCXInKhqv4rosztwIPAylwGGuZxoAhYFtD+4wn62JgcssrEmL0eU9Wl4QkiUh+4D7gGeFtEOqvq4tB6VV0FrMpplGFUdS2wNqj9xxP0sTG5Zbe5TE6ISE0Ruc27dVQiIptE5H0RGRgj/y9FZJaIrBKR7SLynYi8KyJXR+RrLSJjReQrESkVkR+9fTwtIo3TjVtVN6rqtcBzQH3cf9rh+4/aLpAofhFp5d3e6ua9D7/NVhy2naXeUk9Exng/7xCRkd76uLeaRKSZd6tujXd85onI+VHyDfW2MzTGdsrEBdztvf1nePyJjo23bqCIvCciG724/isit4tIzSh5Q8egtog8IiLLvGP6lYjcKiISLWaTW3ZlYrLOa2+YgfviXAw8ARQAA4BJItJBVe8Iyz8MeAZYDbyB+8+7KXAscAnwpJfvIGAO7pbUNOA1YD/gUGAI7hbQugz9GvcCFwF9RKSeqm6K8/v6iX8DcA8wFDjE+zlkacQmawD/ABoBbwObgG98xNwQmO3t669AA2Ag8IKINFfVR3xsI5bHgF/hzumEKDHHJCIP4G6BrQVeBLYAvYEHgF+IyBmq+lNEseq4z1Az4C1gp7f/B3Hn/B5MsFTVFltSWgB1H6GE+W738k4DqoWlN8V9CSnQNSx9HrAdaBplW03Cfr7GK/vbKPlqA7V8/h6hGFolyLfcy3daWNr4yLJ+4/feF8c7hmGxzQRqR1k/0lvfPdq5AV4GqoSlHwr8CPwEtA5LH+rlHxrnXBf72XeCY9PFS1sGHBiWXg1X8SpwR4xjMC38nHqfnw3eUj3ov4d8X+w2l8mFS3FfBjeo6s5QoqquwbVHAFweUWYnsCNyQ+raCCKVRsm3VVXLpKcp1JC8v4+8ycTvx42qujXJMruAW1V1d9j+vwH+iPtPf0iKsaTjUu/1flVdHRbXTuBGYDdlPwsh14afU+/z8zru9mPb7IRr/LLKxGSViNQFDge+07CG6zD/8F5/Fpb2Au422EIReVREfiUi0b7A/467RfKEiLwmIsNE5Ogs3kMPbTfRo7x+4/drG/CfFMot8yqPSMXe68+irMu2473Xf0SuUNUvgBXAod6DD+E2qupXUba33HttmLkQTSqsMjHZFvpSiPVUTyi9QShBVccAFwPfAtcCU4DvReSfItIxLN+3QCdgMtAT106xAPhWRK7N5C/haea9/hAvk9/4k7BGvfs6Sfo+RnroiiDyCzsXkv48eDbEyB+60q2aTlAmfVaZmGzb6L0eGGP9QRH5AFDV51S1M9AYOBt4FjgVmBH+X76qLlLVQV6+jsBtuM/1H0Tkskz9EiJyOHAw7strXqL8fuP3KdVOjQfESA+di/BjHroVVuahHBGJ/GJPR0qfB1P+WWViskpVNwP/A5qLyBFRspzmvX4ao/wGVZ2mqlfgGnQb4b6UI/PtVNV5qvoQUOgl/yrd+MPc5b2+4f1OvviIfxeAiGTjP+uW0R7LBbp7r5+Fpa33XltEyR/ramqX95pM7KF9do9cEVZhf6Oqsa5ETDlllYnJhb/g2hseCf/SFJEmwJ1heULpp8Vo92jqvZZ4+U6Icm8d9v5HXpJu4F7/jj/iGqs34K58EpXxFb8n9Ohyy7QCja4q8JCI7Pk7F5FDcbfedgLPh+Wdi7s6OV9ECsLyNwIejrH9VGIPnecR4Vdo3udiFO476dkktmfKCetnYtImIuPjrL4a9yXRG+gH/FtEpuEaqM/DfcE+rKofhJWZAmwRkY9xj4UKcApwIu4W00wv3xDg1yLyAe7qZz1wGNAX92juY0n+KteJyAZvf6HhVE7FPWb8BXCh10iciN/4AWbhjsNk77iUAt+q6sQkY4/mP8DPgXki8jZ7+5k0AG5R1f+FMqrqKhF5AXdM54vIVNwxOAt4j+iN9f/EVUC/F5H2eFc3qnp/rIBUdbaIPAzcAiwQkVeBrbjPR3vgAyCd/i8mKEE/m2xLxV3Y25ch3tLAy7sfcAeugbwU2Iz74iiMst0rcV/IX+P+i/8Rd3vkFqBuWL6fA08B//bylAJf4TrotU/i91gaEfMOb3v/BSbiOlfWiFF2PGX7UviK38tbFddZ72tvv/v05/BiWxon9pHE7mdSjHto4HlgDe6psE+B82Nsqybui3wFrh/KV7g+QtUi4worcyEw3zv2+/Q7inZswtYN9s7/Zi+uz4H/A/aLcX6iHoNYv78tuV/EOyHGGGNMyqzNxBhjTNqsMjHGGJM2q0yMMcakzSoTY4wxacuLR4ObNGmirVq1CjqMvLdk3RIA2ja2MfkqjE3unFHPzllKlnjHr23FPH7z5s1bq6q+RmzIi8qkVatWzJ07N+gw8l738d0BKB5aHGgcJgkzu7vXnsVBRlFxde/uXouLg4wiZSLyrd+8dpvLGGNM2vLiysSUDyNOHRF0CCZZ7e2cpWVE/hw/q0xMzvRs3TPoEEyyDrRzlpae+XP87DaXyZn5q+czf/X8oMMwyVg/3y0mNfPnuyUP2JWJyZnrpl8HWAN8hTLPnTNrgE/Rdd7xq6AN8MmwKxNjjDFpsysTY0zKNm3axJo1a9ixY0fQoZRPd9/tXhctCjaOCNWrV6dp06bUq1cvY9u0ysQYk5JNmzbx/fff07x5c2rVqkX0+cDyXBXv5k856rSoqpSWlrJy5UqAjFUodpvLGJOSNWvW0Lx5cwoKCqwiqUBEhIKCApo3b86aNWsytl27Mkmg1W1Tgw4hUEsfPDtj23qgxwMZ25bJkeNin7MdO3ZQq1atHAZTATVvHnQEMdWqVSujtyetMjE507VF16BDMMnaP/45syuSBOrUCTqCmDJ97uw2l8mZ2ctnM3v57KDDMMn4YbZbTGq2bHFLHrArE5Mzd8y6A7B+JhXKv905s34mKfIauctTA3y22JWJMSavjRw5EhGJujz//PMUFxfved+gQYM95ZYuXYqI8Oabb8bcdvchQ5B27RARHn/88Vz8OoGxKxNjTN6rX78+06dPL5N++OGHs2DBAgBeeOEF2rRpk9R2n7z7bjZt2UKXwYMzEmd5ZpWJMSbvVatWjc6dO8fNc+yxx9K+ffuktnvU4YenE1aFYre5jDHGpM2uTEzOPHbmY0GHYJJ1Qv6cs507d5ZJq1Ytza/IFi3SK1+BWGVicqbDgR2CDsEkq2Hy5yw0PXO4gUcP5OoTr6ZkRwlnvXBWmfVDOwxlaIehrC1Zy4CXB5RZf1XHqxjUfhDLNy5nyJQhZdbf2OVG+rbtm3SsIevWraN69epl0r/55puUtwlAQUF65SuQnFcmInIU8CegC7AB+DNwj6ru8lG2P3A70B4oAeYA56rq1uxFbDJl5tczAZskq0JZ7c5ZZZ8kq379+sycObNMerNmzVi6dGnqG960KfWyFUxOKxMRaQjMBBYC/YDDgNG4tpu481uKyOXA48DDwM1AQ+B07Oqqwrj/vfsBq0wqlAXunCVTmcTrR1RQvSDu+iYFTeKub1G/RVb6KVWrVo2OHTtmfLusWpX5bZZTuf4ivhKoBfRX1U3AOyJSDxgpIg97aWWISBPgUeAaVR0XtmpK1iM2xhiTUK6f5uoNzIioNIpwFUy3OOUGeq8TshWYMcaY1OX6yqQd8I/wBFVdJiIl3ro3YpT7ObAEuExE/g84APgUuF5VbeAgY0xadu7cyccff1wmvYWPp7E+/PBDtm3btk9aq1atsnPbrBzLdWXSENfoHmm9ty6WA4G2uHaVW4B13ut0ETlCVb+PLCAiw4BhAC1btkwzbGNMZbZx40a6dOlSJv2+++7j5JNPjlv2wQcfLJN28cUXM378+EyFVyFUlMZrAeoA56nqdAARmQ18CwwH7owsoKpjgbEAHTt21NyFamJ5ps8zQYdgktWp8p+zkSNHMnLkyJjri4uLAdi1axe7du2iatWqgLv6UI3/1bLr4IMT5qksct1msh6oHyW9obcuXjkFikMJXrvLPOCoDMZnsqhtk7a0bVL5R0+tVOq1dYuhQ4cONG7cOKkyPc4+m+p162YpovIl11cmi3FtI3uISAugwFsXyyLc1UnkbC4C7M5kgCZ73ljimsTS6VxmcmyF14x5cP6esxNOOIE5c+YAyfeIf+aRR9i8ZQvUrcshhxySjfDKjVxXJm8BN4tIXVXd7KUNAkqBd+OUexO4GzgNmAYgIvWBE4BR2QvXZNLoj0YDVplUKIvdOcvnyqRu3bopN6a3rVcP6tWz+Uyy4GlgOzBZRHp6jeQjgTHhjwuLyFci8mzovarOBV4HnhWRi0XkbODvwA7giVz+AsYYY8rKaWWiquuBHkBV3GPA9+A6I94dkbWalyfchcDfgDHAq7iK5HRvm8YYYwKU86e5VHUhbhiUeHlaRUnbAlzlLcYYY8oRm8/EGGNM2ipKPxNTCUw8Z2LQIZhkdbFzlpZDDw06gpyxysTkTIv6+TNRUKVR285ZWmrUCDqCnLHbXCZnJi2YxKQFk4IOwyTj20luMan58Ue35AGrTEzOPDX3KZ6a+1TQYZhkfPmUWyqhvn37cswxx8RcP3z4cBo0aMD27dvjbmfp0qWIyJ6lTp06HHfccfz5z3+GH35wS5iVK1dSt25d/ve//+1JExEef/zxmPvo06cP9913n8/fLBhWmRhj8lJhYSELFixg4cKFZdbt2rWLV199lf79+1OzZk1f2xs1ahQfffQRU6ZM4bjjjuOKK67g+b//vUy++++/nz59+nDYYYf5jvXWW29lzJgxbNgQbZzc8sEqE2NMXurXrx8FBQW89NJLZdb985//5Pvvv6ewsND39tq2bUvnzp3p1asXEyZM4Mgjj+S5v/1tnzybNm1iwoQJXHrppUnFesopp9C4cWMmTiy/D0RYZWKMyUu1a9emb9++TJpUtk2oqKiIpk2b0qxZMwYPHkyLFi0oKCjg6KOP5rHHHmP37vhDAooIxxxzDMtXr94n/eWXX6ZWrVqcfnrcrnZRnXvuuTz33HNJl8sVe5rLGJNZM7uXTWs5ENpcDTtLoPissutbD3XLtrXwwYCy64+4Cg4ZBFuXw0dDyq5vd2NK44cVFhYyadIk5s2bxwknnADAjh07mDx5MhdccAGrV6+mbdu2XHDBBdStW5f58+dz9913U1payu233x5328uWLePQgw/eJ23WrFl06tRpzzD2yejatSuPPPII69evp2HDeNM/BcMqE5Mzrw58NegQTLJOrtznrHfv3jRo0ICioqI9lcmMGTNYv349hYWFdO3alR49egCgqpx88smUlJQwbty4MpXJ7t272blzJ5s3b2bChAl8+umnvPPWW9C69Z488+bNo1+/finFetxxx6GqzJ07l169eqX4G2ePVSYmZ5oUNAk6BJOs/VI4Zz2LY6+rVhB//X5N4q+v3SL++iTVqFGD/v378/LLL/Pwww8jIkyaNIlDDjmELl26sG3bNn7/+9/zwgsvsGzZMnbs2LGn7M6dO/cZkj6ykvjDH/7AqRG3s1avXk2TJqn9HYTKrY64dVZeWJuJyZnx88czfv74oMMwyfh6vFsqscLCQpYtW8ZHH33Etm3beP311xk8eDAiwq233sqoUaMYNmwY06ZNY86cOYwYMQKgzLzvjz76KHPmzGHq1Kl07dqVm266iX8XF8PatXvybNu2zffTYZFC5SL3W17YlYnJmVBFMrTD0EDjMEkIVSSthwYZRVaddtppHHDAARQVFbFq1So2b9685ymuV155hWuuuYZbbrllT/6pU6dG3c7hhx++Z96TLl26cMQRR3DbXXfx1rhx4F1VNGrUKOXHe0PlGjVqlFL5bLMrE2NMXqtatSoDBw7klVde4cUXX+TII4/kuOOOA6C0tHSfK4ldu3ZRVFSUcJsNGzbk1ltvZfr77/OfJUv2pLdt25ZvvvkmpTiXLl0KQJs2bVIqn21WmRhj8l5hYSGrV69mypQp+/Qt6dWrF0888QQTJ05k6tSp9O3bN2GP+JCrrrqKRvXr88ize+b546STTmLevHlR88+fP59XX311n+Xdd/dOQDt37lzq16/P0UcfneJvmV1WmRhj8l6XLl1o1aoVqrpPZfKnP/2JU045hd/85jdceumltG/fPuEjwSF16tThtxddRNG0aSxfvhyA/v37s3DhQpYtW1Ym/7PPPst55523z3L33XvnDZw+fTrnnHMOVaqUz69tazMxxhiIevvpgAMOYMqUKWXSr7jiij0/hyqhaO76zW+46ze/gRZu9OWOHTtyzDHHMGnSJG6++eY9+WKVD9m4cSMzZsxg5syZvn6XIFhlYnJm2gXTgg7BJKu7nbO0HH54maQRI0Zw8803c/311+/zaHE8Tz31FJ07d+bkk0/OdIQZY5WJyZmC6gVBh2CSVc3OWVqi9HQfMGAAX3/9NStXruSQQw7xtZn69evzxz/+MdPRZVTOb76JyFEiMktESkTkOxG5V0Tiji0gIq1ERKMsiR+rMOXGk3Oe5Mk5TwYdhknGF0+6xaRmzRq3hAn1X/FbkYBrzD/22GMzHV1G5fTKREQaAjOBhUA/4DBgNK5SG+FjEzcBH4a9Xxsroyl/Xv78ZQCuPvHqgCMxvi1z54w2ds5Ssn69e23aNNg4ciDXt7muBGoB/VV1E/COiNQDRorIw15aPEtU9eOsR2mMMSYpSd3mEpEaItJcRA7zrjKS1RuYEVFpFOEqmG4pbM8YY0w5kLAyEZGjReQhEZkHbAGWAV8Aa0VkjYj8TUQuFJFaPvbXDlgcnqCqy4ASb10ifxWRXSKySkTG+NynMcaYLIt5m0tETgLuB04F5gDvAn/EtVNsBxoArYCOwKPAn0RkDPCoqm6JsdmGQLSBadZ762LZDjwBvA1sAroDt+LaXKKO5ywiw4BhAC1btoyzaWOMMemK12YyGVd5DFHVFfE24j2N1RO4zku6LzPhOaq6ChgellQsIt8DT4rIcar67yhlxgJjATp27Bi/R5DJieKhxUGHYJKVweHe81LbtkFHkDPxbnMdoqq/S1SRAKjqLlWdoaq9gUfiZF0P1I+S3tBbl4zQrD0nJFnOGGPKeO211zj99NNp0KABNWvWpE2bNtxwww18+OGHiAivvfZa1HLff/891apV46GHHvK1n6FDhyIiiAhVqlTh4IMPprCwcM9AjpGuvfZaLrnkkj3vR44cGXdOlLlz59KoUSM2btzoK55MiVmZqGpKg+YnKLeYiLYREWkBFBDRluJnVxGvppwbNXsUo2aPCjoMk4xFo9xSyd14440MHDiQ1q1bM3HiRN5++22uv/56Zs2axahRozjiiCNijhb8yiuvsHv3bgYPHlx25erVbonQrl07PvroIz744APuvfdeiouLOeuss/jpp5/2ybd8+XLGjRvHrbfe6vt36dixIz/72c949NFHfZfJhKQfDRaRari2iNMAAf4JPKOqO30Ufwu4WUTqqupmL20QUIprk0lGaKLo6ENwmnLnzS/eBOCmrjcFHInxbaU7ZxxZec/ZG2+8wZgxY3j22We59NJL96R369aNYcOG8fbbb/Pxxx/zyCOPsGXLFurUqbNP+aKiIrp06RK9E2Lo6uDAA/dJrl27Np07dwbc3O4FBQUUFhYyd+5cunbtuiff008/zfHHH0+7dn6eT9rrkksu4aabbmLEiBG+h2xJVyo94P8AXAQUA/8CbsE1jvvxNK4xfbKI9PQayUcCY8IfFxaRr0Tk2bD3I0VktIj098rdi2v0n6yq/0nhdzDGGMDNkHj88cfvU5GEVK1ald69e1NYWEhpaSmvv/76PuuXL1/O7Nmz94w0PHr0aE488UTq16/PAQccQN8rr+Srb79NGENo/pTQ6MIhzz33HAMGDIhWJK5f/vKX/Pjjj8yYMSPpsqmKWZmIyPExVg0AzlDVJ1T1YeBq4Dw/O1PV9UAPoCrwBnAPrlK4OyJrNS9PyGJcP5S/AtOA83FtM+f72a8xxkSzY8cOZs+ezZlnnhk3X7t27ejQoUOZW12TJk2iSpUqnHee+wpcsWIFw4cP5/XXX2fcuHHs2r2broWFCdsvQkPSH3rooXvSlixZwooVK/a5UvGrXr16HH300TkdZTje9c9bIvJ34A5V/SEsfTXQC3hNRAR3u2uV3x2q6kLg9AR5WkW8L8J1bjTGlHfdu5dNGzgQrr4aSkrgrLPKrh861C1r10K0/8SvugoGDYLly2HIkLLrb7wR+vZNOtR169axfft2X90HCgsLufPOO1m/fj0NG7qeDEVFRZx++ukccMABAPu0U+zatYtehxxC065def3117nooov22d7OnTtRVRYtWsRtt93GmWeeSadOnfasD02i1b59+6R/L3BXO//6179SKpuKeLe52uI6Ey4SkZtFpLqXPhzXp2QN7gms84GrshumqQxqVa9FrerWz7RCqVrLLZWc+784vsGDB7Njx44985v873//Y968eftMpvXxxx/Tq1cvGjduTLVq1Sjo0IEtJSV88cUX+2xr3rx5VK9enRo1anDcccexadMmXnrppX3yrF69mv3224/atWun9Ds1adKE1VEa/7Ml5pWJqm4AfisiT+MGYxwmIjeq6t9FpBV7n8paoqr+5rE0ee2tC94KOgSTrNNSOGfFxbHXFRTEX9+kSfz1LVrEX5+kxo0bU7NmzagzH0Zq2bIlXbt2paioiEsvvZSioiJq1qxJ//79AXer6owzzqBTp04888wzNGvWjBo1anD22Wezbdu+D7keeeSRPPfcc+zYsYMPPviAO+64g1//+tdMmjRpT55t27btM/98smrWrFlmv9mUsJlfVRcBZ4nI2cAoEbkG+K01fBtjKrrq1atz0kknMWPGDO6///6E+QsLC/ntb3/LDz/8QFFREb1796Z+fdd1bvr06ZSUlPD666/vuZrYuXMnP/74Y5ntFBQU0LFjR8BNGbxt2zbuuusubrjhBn7+858D0KhRIzZt2sTu3btTmqp3w4YNNGrUKOlyqfIzNlctEamvqlOBY3CP974rIo+LSO4iNRXefe/ex33vZnRwBJNt/73PLZXYddddx9y5c5kwYUKZdbt372b69Ol73oca2u+55x4WLFiwzy2u0tJSqlSpss+juC8/8ww7dybuNXHjjTfSpEmTfTo+tm3bFlXlWx9Pg0WzdOlS2rRpk1LZVMQbm+twYALQBVAR+Rq4UlXHiMhE3JApi0XkPuBJVd2Vk4hNhTXrm1kA3NntzoAjMb59784Zx1Tec9a3b19uuOEGLrvsMj788EP69etHnTp1WLx4MU8//TStWrXa87RX06ZN6dGjB08++SR16tShb1ij/+mnn86uXbu45JJLuOyyy/j8888Z9eCDNKhXL2EMBQUFXH/99dx55518+eWXHHHEEXTq1Ilq1aoxb968fZ7yAvjpp5949dVXy2ynW7du7L///oDrCZ9MZ8d0xbsyeQ5YChyIG9RxAq5/yH6q+oOqXol7qusc4L/ZDtQYY7Jl9OjRTJo0iS+//JLzzz+fXr16MXr0aHr06MFTTz21T97CwkJUlX79+lGr1t6HE4455hjGjx/PJ598Qp8+fXjxxRd55bHHqB/RyTGW4cOHU69ePUaNciMO1K5dm1/84he89VbZdqvNmzdz3nnnlVk+//xzAD777DN++OGHPe05uSCq0UcjEZENwHmq+o73vhFuxOA2qvpVRN5zVHVKtoNNVceOHXXu3LkplW1129QMR1OxLH3w7Ixtq/v47oAN+FihzOzuXqMM+Lho0SKOPPLInIZT4SxZ4l5THPBxypQpXH755Xz33XdJNcbffvvtzJkzJ2E/k0TnUETmqWpHP/uMd2XyNvCQiJwrImcBfwH+5y37KM8ViTHGVFT9+vWjWbNmTJw40XeZrVu3Mm7cOEaM8DMTeubEe5rrUty87LcBNXBjYPXUWJcyxiTQuKBx0CGYZNW0c5aWqlUT54mjSpUqjBs3jiWhKxwfli1bxl133UX3aJ1HsyheP5MtuIrEmIx4bWD0IbxNOXaKnbO0HH542pvo3LnznkEh/TjyyCMDuf0Yb2yuxF1CM1jOGGNMxRWvzeQLEblcRHz15ReRE0TkOexqxsRw+8zbuX3m7UGHYZIx/3a3xGB3vRNYscIt5VCmz128NpNbcaP6/kFE3gZmAwvYdw74Q3EzHZ4JtAD+jGuoN6aMj1Z8FHQIJllrY5+z6tWrU1paSkFBQQ4DqmC2bg06gphKS0upXr164ow+xWszmSwiU3Bzu1+Em9/9IPbObCjAT7iG+WeAiaq6JmORGWPKtaZNm7Jy5UqaN29OrVq1fA2WaIKnqpSWlrJy5co9ox1nQtyxubwnt97xFkTkIFwnxv2AH4GlNsijMfmpntez+7vvvmPHjh0BR1NOhUbt3b072DgiVK9enQMOOGDPOcyEpOZzVNVVJDF3iTGmcqtXr15Gv5Aqnau82TkyONJxeZWbyYGNAQ6ud3DQIZhkFdg5S8vB+XP8rDIxOfN8/+eDDsEkq6uds7Q8nz/HL/lB8o0xxpgIOa9MROQoEZklIiUi8p2I3CsivsccEJEqIjJXRFRE+mQzVpNZ102/juumXxd0GCYZ865zi0nNdde5JQ/4us0lIseoatrDzItIQ2AmsBDoBxyGmxK4Cm4cMD8uB9S0uIgAAB6dSURBVPLnRmQlMn/1/KBDMMlab+csLfPz5/j5vTL5t4jMEZGrRKRBGvu7EqgF9FfVd1T1aVzHyBtEJOEjIV5l9Dvg/9KIwRhjTIb5rUxOx11NPAx8JyIviUivFMbh6g3MUNVNYWlFuAqmm4/y9wEfArOS3K8xxpgs8lWZqGqxql6M67A4HGgOzAC+FZH7ROQwn/trByyO2PYyoMRbF5OIHIsbFv8mn/syxhiTI0k1wKvqVlX9i6qeCrTFTet7B25QyHdF5JwEm2gIbIiSvt5bF8+fgMcjZ3mMRUSGeQ31c3/44Qc/RUyWtWnchjaN2wQdhklG3TZuMalp08YteSDpfiYi0goYihuvqwUwDfgb8Atgkog8oarXZy5EEJHBuMqrr98yqjoWGAtu2t5MxmNSM7bv2KBDMMn6uZ2ztIzNn+Pn68pERApE5CIR+SfwFXABMA5oqap9VfVZVR0I/Bq4LM6m1gP1o6Q39NZF23d14BHgIaCK9wBAqLG+tojU9fM7GGOMyR6/Vybf4yqeybipe4tj5JsDrIuzncVEtI2ISAuggIi2lDC1cY8Cj/GWcEW4OenTn87MZN2wN4YBdoVSoXzizpldoaRomHf88uAKxW9lcgvwoqpujJdJVRfg5jiJ5S3gZhGpq6qbvbRBQCnwbowyW4DTItIOBF7Ctdf8I0Hsppz4Yt0XQYdgkrXZzllavsif4+e3AX5/3BVCGSJykIjc5XM7T+Mm1posIj1FZBgwEhgT/riwiHwlIs8CqOpO72myPQvwsZf1v6r6ic99G2OMyRK/lcndxO513sxbn5Cqrgd6AFWBN3AdFh+NUr6al8cYY0wF4Pc2l7B3hsVIBxOj8TwaVV2I6wQZL0+rBOuXejEZY4wpB2JWJiJyMXCx91aBp0RkU0S2/YBjgLezE56pTDoc2CHoEEyyGto5S0uH/Dl+8a5MStj7ZJYAG3FT9Yb7Cdeo/mTmQzOVzWNnPhZ0CCZZJ9g5S8tj+XP8YlYmqvoK8AqAiPwVuE9Vv85VYMYYYyoOX20mqnpJtgMxld+Fky8EbMbFCmW2O2c242KKLvSOXx7MuBivzeRh4I+qusL7OR5V1VszG5qpbFZsWhF0CCZZJXbO0rIif45fvCuT84AXgBXez/EoYJWJMcbkqXhtJodG+9kYY4yJlPM54I0xxlQ+8dpMzkpmQ6o6Lf1wTGXW5eAuQYdgktXEzllauuTP8YvXZvImri3ET09zxYY/MQn8vufvgw7BJKuDnbO0/D5/jl+8ysTaSYwxxvgSrwH+21wGYiq/c18+F4DXBr4WcCTGt/fdOeMUO2cpOdc7fq9V/uMXr82kQFVLQj8n2lAorzGxrCuJN2+aKZe22zlLy7r8OX7xbnNtFpEuqvov3ARVieZRtzYTY4zJU/Eqk0txU+KGfk5UmRhjjMlT8dpMJoT9PD4n0RhjjKmQ/E6OBYCINADaAwcBq4AFqrohG4GZyqfHoT2CDsEk6wA7Z2npkT/Hz1dlIiLVgN8BvwHCG+NLRORJ4P9UdUcW4jOVyJ3d7gw6BJOsY+ycpeXO/Dl+fq9MxgDDgHuBycAaoClwLjACN+PitdkI0BhjTPnnd2yuIcAdqvqAqi5W1R+919/hKpMhfncoIkeJyCwRKRGR70TkXhGJ+ySYiBwtItO9/NtFZJmI/FlEDvK7XxO83i/0pvcLvYMOwyTjn73dYlLTu7db8oDfK5PdwOcx1i3A55NeItIQmAksBPoBhwGjcZXaiDhF6wPfAM8B3+F6598NnCAiJ6rqTj/7N8Eq3VEadAgmWbvsnKWlNH+On9/KZCJwOTAjyrorAL/TiF0J1AL6q+om4B0RqQeMFJGHvbQyVHU2MDssqVhEVgBvA8cCn/rcvzHGmCyI1wP+6rC3S4EBIvI58Hf2tpn0A+oCo3zurzcwI6LSKAIeAroBb/iOHEJdS2skUcYYY0wWxLsyeTxKWjPgyCjpY4A/+NhfO+Af4QmqukxESrx1cSsTEamCi/lQ4EFgDvAvH/s1xhiTRfE6LWZj4qyGQLR+Keu9dYlMA37h/TwPOEtVd0fLKCLDcE+g0bJly+QjNRnXp02foEMwyWpu5ywtffLn+CXVabEcuAZoBByBa7B/S0ROUtVtkRlVdSwwFqBjx442FEw5cFPXm4IOwSTrSDtnabkpf45fsj3gDwba4PqV7MPnTIvrcU9mRWrorYtLVb/0fvxERN7HPeF1PvAXH/s2xhiTJX57wNcFXgbOCCV5r+H/8fsZNXgxrm0kfNstcL3qF/uJJURVvxWRH4HWyZQzwek+vjsAxUOLA43DJGFmd/faszjIKCqu7t3da3FxkFHkhN92kd8DLYFTcBXJOUB34Fnc1UFnn9t5C/iFVzmFDAJKgXd9bgMAEWkLNPb2b4wxJkB+b3OdhWuj+MR7/52qzgHeE5HRwM3AQB/beRo37MpkEXkId1UxEhgT/riwiHwFvKuql3nvRwE7vf1vwD1RdgtuiPwin7+DMcaYLPFbmRwALFfVXSKyFdcIHjIN8DUnpaquF5EeuMeO38BVDI/iKpTIuMJvm83FNb4Pw7XXLPP2+XtV3erzdzDGGJMlfiuT5UAT7+cvgT7s7Q3/c6DM01SxqOpC4PQEeVpFvC/CrkCMMabc8luZvAP0BKbgriQmiMgJwHbgVNz4WsbENfBoP3dCTbnS0s5ZWgbmz/HzW5ncijePiapOFJEtwADcOFvDgWeyE56pTK4+8erEmUz50sbOWVquzp/j56syUdUSoCTs/RTcVYoxvpXscB+hguoFCXKacmOn92dfzc5ZSkq841dQ+Y9fsp0W2wIn4qbt/Q6Yq6pLshGYqXzOeuEswPqZVCjF7pxZP5MUneUdvzzoZ+K302I9YBxuZsUqwBagDrBbRCYDl8caPt4YY0zl57fT4pO43u8XAbVVtR5QG7gY6OWtN8YYk6f83ubqB1yvqi+GElS1FHhBRApwQ9AbY4zJU36vTLYAq2Ks+w6wjoPGGJPH/F6ZPAHcJCL/8K5IAPCuSm7CbnMZH4Z2GBp0CCZZrYcGHUHFNnRo0BHkTLxpex+OSDoCWC4i77B32t5euEEa52YtQlNpWGVSAVllkh6rTAA4L+L9Dm8JHyF4s/d6Lm6wR2NiWluyFoAmBU0S5DTlxjZ3ztjPzllK1nrHr0nlP37xpu09NJeBmMpvwMsDAOtnUqF84M6Z9TNJ0QDv+OVBP5NszPNujDEmz/iuTESktYg8JSL/FZGV3uuTImIzHRpjTJ7z2wP+BOCfuKHm3wS+x81xci5wgYicpqqfZi1KY4wx5ZrfR4NHAZ8Bvb1BH4E9jwZP89bHnaPEGGNM5eW3MukEDAyvSMCNJuxNqTsp45GZSueqjlcFHYJJ1hF2ztJyVf4cP7+VSSnQOMa6RiQx06LJX4PaDwo6BJOsQ+ycpWVQ/hw/vw3wU4EHReTk8ETv/e9x87kbE9fyjctZvnF50GGYZGxd7haTmuXL3ZIH/FYmNwBfA++KyCoR+beIrALeBb4BbvS7QxE5SkRmiUiJiHwnIveKSNUEZU4Ukb+KyFdeuSUicreI7Od3vyZ4Q6YMYciUIUGHYZLx0RC3mNQMGeKWPOB3psV1wMkiciZ7J8daBXyiqm/73ZmINARmAgtxIxEfhps/vgowIk7RQV7eh4AvgWOB+7zXc/3u3xhjTHYkrEy8//7/A1yrqtOB6Wns70rcvPH9vcm03vEm3hopIg/HmWDrQVVdG/a+WES2Ac+IyCGq+m0aMRljjElTwspEVbeJSANgdwb21xuYEVFpFOGuOLoRo+0loiIJ+cx7bQZYZWIqpVa3TQ10/0Wt1wEwOMA4lj54dmD7Nv75bTN5AbgkA/trBywOT1DVZUCJty4ZXXAV3P8yEJcxxpg0+H00eBkwUETmAG/hesBr2HpV1ad8bKchsCFK+npvnS8iciCujWWiqq6JkWcYMAygZcuWfjdtsujGLr6f0zDlxLgfzgk6hIrtxvz5zPutTEZ7rwcBJ0RZr4CfyiRtIlIDeBk3++P1sfKp6lhgLEDHjh01Vj6TO33b9g06BJOkWZt/HnQIFVvf/PnM+32aK1OjC68H6kdJb+iti0tEBHgOOBo4SVUTljHlx5K1SwBo26RtwJEYv1rXXAHA19sPDjiSCmqJ+8zTtvJ/5v1emWTKYiLaRkSkBVBARFtKDI/hHinupap+8pty5Ndv/hqw+UwqkgeaPw7A4K8fDDiSCurX7jOfD/OZ+K5MvNtLQ3HjdO3pZwJMUNWffG7mLeBmEamrqqFZGgfhhmt5N8H+bweG48YI+8Bv3MYYY7LP1+0rETkS11nwCaA9sMt7fQL4SkSO8rm/p4HtwGQR6ek1ko8ExoQ/Luz1dH827P35wAO4W1wrRaRz2LK/z30bY4zJEr9XJmOBjcAp3qO8AIhIS9z8Jk8DpybaiKquF5EewOO4PiUbgEdxFUpkXOFDrJzhvQ71lnCXAON9/RbGGGOywm9l0hEoDK9IwPUREZG7gRf97lBVF5Jg7hNVbRXxfihlKxFjjDHlhN/KZCkQa1DF/XD9UIyJa8Sp8YZfM+XRn9YMDjqEim1E/nzm/VYmtwGjReQbVf0klCginXEDLt6UjeBM5dKzdc+gQzBJ+nBLh6BDqNh65s9n3m9lMgKoB8wWkTXAGqCpt6wD7hCRO0KZVbVTpgM1Fd/81fMB6HCgfUFVFEft9zUAC7e1DjiSCmq++8zTofJ/5v1WJgu8xZiUXTf9OsD6mVQkdzUbC1g/k5Rd5z7z1s/Eo6qZGOTRGGNMJZWpYVKMMcbkMatMjDHGpM0qE2OMMWnL9UCPJo890OOBoEMwSXp49cVBh1CxPZA/n3mrTEzOdG3RNegQTJI+LTky6BAqtq7585m321wmZ2Yvn83s5bODDsMk4fiCRRxfsCjoMCqu2bPdkgfsysTkzB2zXL9W62dScdxy4ATA+pmk7A6vL3ce9DOxKxNjjDFps8rEGGNM2qwyMcYYkzarTIwxxqTNGuBNzjx25mNBh2CSdO93w4IOoWJ7LH8+81aZmJyxoecrHht6Pk15MPR8iN3mMjkz8+uZzPx6ZtBhmCScVGc+J9WZH3QYFdfMmW7JAzm/MhGRo4A/AV2ADcCfgXtUdVecMjWA3wGdcfPR76eqkoNwTQbd/979gM24WJFc07QIsBkXU3a/+8znw4yLOb0yEZGGwExAgX7AvcCNwD0JihYAlwMlQH50JzXGmAok11cmVwK1gP6qugl4R0TqASNF5GEvrQxV3SAijVRVRWQ4cHoOYzbGGJNArttMegMzIiqNIlwF0y1eQVXVbAZmjDEmdbm+MmkH/CM8QVWXiUiJt+6NHMdjsqzVbVP3/Ly6xroyaZXd0gfPDjoEY3Ii15VJQ1yje6T13rqMEZFhwDCAli1bZnLTJkWNdwwPOgSTpDtW2jlLyzPPBB1BzlTafiaqOhYYC9CxY0e7RVYOVNeDgw7BJOnr7XbO0tK2bdAR5Eyu20zWA/WjpDf01plKrKTKJ5RU+SToMEwSetT9hB517Zyl7I033JIHcn1lshjXNrKHiLTAPfq7OMexmBzbVG0KAAU//TzgSIxfV+zvztmszXbOUjJ6tHvt2zfYOHIg11cmbwG/EJG6YWmDgFLg3RzHYowxJkNyXZk8DWwHJotIT6+RfCQwJvxxYRH5SkSeDS8oIr1FZADQwXs/wFsOyV34xhhjosnpbS5VXS8iPYDHcY8BbwAexVUokXFVjUh7CgivOF7xXi8Bxmc6VmOMMf7l/GkuVV1Igh7sqtrKT5oxxpjyodI+GmzKnyY7bgw6BJOk65fbOUvLxIlBR5AzVpmYnKmm+wcdgknSqh12ztLSokXQEeSMzWdicmZr1ffYWvW9oMMwSehT/z361LdzlrJJk9ySB+zKxOTM5qrTAKi969SAIzF+XdjYnbM3N9o5S8lTT7nXQYOCjSMH7MrEGGNM2qwyMcYYkzarTIwxxqTNKhNjjDFpswZ4kzP7/3R70CGYJF31rZ2ztLz6atAR5IxVJiZnqkadfcCUZ+t32TlLS5MmQUeQM3aby+TMlqoz2VJ1ZtBhmCQMaDiTAQ3tnKVs/Hi35AG7MjE5E6pI6uzqGXAkxq9QRfLq+op7zlrdNjWwfRe96OYzGbw4uJEElj54dk72Y1cmxhhj0maViTHGmLRZZWKMMSZtVpkYY4xJmzXAm5xp+tPIoEMwSRr6zcigQ6jQhp43MugQcsYqE5MzVdgv6BBMkrapnbN0bKueP8fPbnOZnNlcdSqbqwb3mKZJ3oWNp3JhYztnqbrw06lc+Gl+HL+cVyYicpSIzBKREhH5TkTuFZGqPsrVF5G/ish6EdkoIi+ISONcxGwyY2vV99la9f2gwzBJ6FP/ffrUt3OWqj6L36fP4vw4fjm9zSUiDYGZwEKgH3AYMBpXqY1IUPxloA1wObAbeAj4G3BKtuI1xhjjT67bTK4EagH9VXUT8I6I1ANGisjDXloZItIFOAPopqrveWkrgU9EpKeq2ngPxhgToFzf5uoNzIioNIpwFUy3BOW+D1UkAKr6L+Abb50xxpgA5boyaQcsDk9Q1WVAibfOdznPogTljDHG5ICoau52JrIDuFlVH4tIXwE8p6p3xCj3DrBVVX8Vkf480FpVu0YpMwwY5r1tCyyJEVYTYG1Sv0huWXzpsfjSV95jtPjSEy++Q1TV1yiVlbafiaqOBcYmyicic1W1Yw5CSonFlx6LL33lPUaLLz2Zii/Xt7nWQ9QZkhp66zJdzhhjTA7kujJZTEQbh4i0AAqI3iYSs5wnVluKMcaYHMp1ZfIW8AsRqRuWNggoBd5NUO5AETk5lCAiHYHW3rp0JLwVFjCLLz0WX/rKe4wWX3oyEl+uG+Ab4josLsB1OmwNjAEeU9URYfm+At5V1cvC0mYARwA3sbfT4hpVtU6LxhgTsJxemajqeqAHUBV4A7gHeBS4OyJrNS9PuEG4q5e/AM8B84BzshmvMcYYf3J6ZWKMMaZyystRg0XkChH5UkS2icg8Eenho8xIEdEoy5kpxlDuB7xMJUYRaRXjOBVlOLbDReQZEfmPiOwSkWKf5XJy/FKJL1fHztvXeSLydxFZKSJbvL+DQh/laorIaBFZIyJbRWSqiLQqR/FFO34fZyG+ASIyW0TWed8jS0RkhIjUSFAuV5+/pONL9/NXafuZxOJ9IJ8GRgIfAJcAb4rIiaq6IEHxjUBk5bEohRjK/YCXacYIrm3rw7D3me60dTRwFvAxUD2JcrkaMDTV+CD7xw7gBtxwRNd72z8LeFFEmqjqn+KU+yMwwCv3A+7v6B0ROUZVt5WD+MB9Tl8Ne785g3GFNAb+ATwCbAA64Y7FgcDwOOVy9flLNT5I9fOnqnm14HrC/yXsfRXgv8DzCcqNBNZmKIbbcf1j6oWl3YIbVqZenHJdAAVODUvr5KX1zPBxSjXGVl48fbJ8HquE/fwqUOyjTC6PXyrx5eTYeftqEiXtReCbOGUOBnYCF4WlNQd+Ai4POj4vjwLDs338Yuz7d7gvbomxPmefvxTjS+vzl1e3uUSkNe6/gpdDaaq6G3iF3A4YWREGvEw1xpzwzluycnb8UowvZ1Q12n+bnwHN4hQ7w3udHLadlbgr/Ewfv1TiC9o6IN5trqAHrE0UX1ryqjJhb8fHyI6Oi4BGIpJoDJoGIrJWRHaIyGci0j+NOMr7gJepxhjyV6+tYJWIjBGRWhmOLxUVZcDQoI5dF+CLOOvbAStUdUtEeq6OX6L4QkaKyE7vb/UvItIoWwGJSFURKRDXB+5a4Cn1/s2PIuefvyTjC0np85dvbSYNvdcNEenrw9b/EKPsV7jbPJ8BdYFfA6+JyLmqOjlGmXhxRMYQiqNhlHQ/5VonGUMiqca4HXgCeBvYBHQHbsW1ufTLbIhJy+XxS0Vgx07cQyi/Ai6Nky3Vz0TafMYHMAHX7eAHoCNwJ3CciHRS1V1ZCG0rUNP7+Tng5jh5g/j8JRNfWp+/Cl+ZiEh94KBE+VQ1rWFXVPX5iP2+AcwG7iLssj/fqeoq9m3gKxaR74EnReQ4Vf13QKGVe0EdO+9prBeB11V1fDb2kY5k4lPVoWFv3xORRcA0oC+uoTvTuuKGg+qE+y54HLg6C/tJle/40v38VYbbXOfhLhMTLbD3CiRy0MiGEesT8i4VJwPHio9HeiNUhAEvM7mv0JM1J6QVUfoq4oChWT123i2gt4BvgQsSZM/58UsyvmimA1uA4zMZV4iqfqqqH6jqGNxtpKtE5LAY2XN+/JKMLxrfn78KX5mo6p9VVRItXvbQ1Unk/cl2wI+qGusWV8zde0uyKsKAl6nGGI1GvAalIg4YmrVjJyIFwJu4Rtk+qlqSoMhioIWI1I5Iz8rxSyG+MsLaB3Lx2fvUez00xvqgP3+J4ovG9/Gr8JVJMlT1a1wD3nmhNBGp4r1PasBIERHgXODfKdyLLY8DXmYqxmgGeK/zMhFYGnJ5/DIlK8dORKrhnmI8AjhTVdf4KPa297pnGCMRaYbrI5HR45difNG2cyZQh9x89k7yXr+JsT7oz1+i+KLx//nL9rPN5W0BCoFduI53pwHjcV+Q7cPydMM9T98tLO1d3GXiGbg/pmm4Tke/TCGGhsAq4B2gJ25GyC3A/RH5vgKejUibAXwN9Mc1SC4B3s/CcUopRlx/nNFefD2Be73j+1qG4yvwPugDgI+Az8PeF5SD45d0fLk6dt6+xuL+27wW6Byx1PTyzAJmRZR7BteJbQiuA+/HwJfAfkHH531GxwIDgdNxne82AJ8AVTMc33Rv+72974R7vL+Polh/Gzn+/CUdX7qfv4z+AhVlAa7wDuR23KVfj4j13b0PcvewtGe9D0Ep7gmJ94HeacRwFK6HainuS/u+yA88sBQYH5HWAPir90eyCdcwWaaDV4aOU9IxAoOBubjRAn7yjvO9oS+ADMbWir23GSOXVkEfv1Tiy9WxC9t3oviKiehsiXsyaAzuaamtuH+qDi0P8eEGkf0Q159iB7Ac12O/fhbiuw83+vkW77P0KXANUD3W30aOP39Jx5fu588GejTGGJO2vGozMcYYkx1WmRhjjEmbVSbGGGPSZpWJMcaYtFllYowxJm1WmRhjjEmbVSYm74ibgjkbsxemRUTGi8jcFMqJiMwXkYtjrC+OkT7Am8412bHljCnDKhNjKr6BQCNcB7hkTAYE15vdmLRYZWJMxXctMFFVd4QSRORAEZkkIuuAbiKyWURmi0j7UB51s0E+h+sZbUxarDIxJgoROV1EPhGRbSLyvYg8KSJ1wtZXF5FRIrJMRLaLyHciMkVEanjrG4jIn730bV6+cVmI83DcnBWvRqwahxtjbjhukL7BuOl1Iyexeg04XkSOznRsJr9U+MmxjMk074t1Om6Qy3OBFsCDuNFdz/Sy3Y6bX+M23CisBwJnAaH2hzG4L/nrgdXeNk7NQrg9cGNkRU5c1A14SFVfEpFfq+pUYGpkYVVdJCLrcQP7fZ6F+EyesMrEmLLuxE3G9Ev1phcQkR+BSSLSRVU/ws1c96KqTggr93LYz52AJ1R1UljaPrN1ZsgJwCLvllW4VUAXEakZpUyk/+DiNSZldpvLmLI6AVN033lqXsNNSxCai2I+MFREbhGRY735bcLNB24WkatFpE0WYz0QNyR8pGtx81esBNqLyG1xZthb623HmJRZZWJMWQcB34cneBXLOtxTUwD3A0/g5tP+N7BcRH4bVmQ4bs7xu4AlIvKliAzOQqz74aZS2IeqzsDdlrvWW38p8LmInBeZ11u/XxZiM3nEKhNjyloFNA1P8PpiNAZ+BFDVbap6l6q2AtoAk4DHvJn9UNUNqnqtqh4IHIeboOkFETkqw7H+iJsjowxVXa+qL+Imr2qHazN5MErWBt52jEmZVSbGlPUJcE5EZ77+uDbGDyIzq+qXuFnttuMmFItc/x/gZtzfW7Q5wNOxhChzekfedvPaVD7FVYiRWuGmszYmZdYAb/JVDREZECX9XdwtrM+Av4nIU8DBwEPADK/xHRGZgnvk9jPcTJQDcH9P73nrPwCm4Ga7U9zsnluBfyWIq2GMuKapakmU9A+Bu0Rkf1X9ISz9MxF51Iuvpoj0BX4DzAwvLCK1cRXcnQniMiYuq0xMvqoLvBIl/TRVLRaR3sADuF7im4CXgFvC8s0GBrH3imMhcK6qhoZD+QgYivuvfxfuS723qq5IEFfrGHEdiptmNVIx7hbVmcDEsPQiXHtOG9xtrGdxjzpfF1H+DKAENze5MSmzaXuNqeBE5A/A4ap6doz1xaraPca6l4Ctqnp5FkM0ecCuTIyp+B4BvhCRNqrqu+1DRFoA/YBjsxaZyRvWAG9MBefdOrsU90hzNONjpB8MXKmqX2UjLpNf7DaXMcaYtNmViTHGmLRZZWKMMSZtVpkYY4xJm1Umxhhj0maViTHGmLT9PyHeEieaz5BcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot loss PDF, expected loss, var, and cvar\n",
    "plt.bar(losses, pdf)\n",
    "plt.axvline(expected_loss, color='green', linestyle='--', label='E[L]')\n",
    "plt.axvline(exact_var, color='orange', linestyle='--', label='VaR(L)')\n",
    "plt.axvline(exact_cvar, color='red', linestyle='--', label='CVaR(L)')\n",
    "plt.legend(fontsize=15)\n",
    "plt.xlabel('Loss L ($)', size=15)\n",
    "plt.ylabel('probability (%)', size=15)\n",
    "plt.title('Loss Distribution', size=20)\n",
    "plt.xticks(size=15)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAElCAYAAADZb/T+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3wc1bXA8d+RZEmW3OQmd8vdmBLANtjg3gikUAI4AUIHEyCATUjhpRB4LyQk2BCqSeiQGEgMJDR32eAC2HSMu+XeLRf1dt4fM5LW69VqV9rdWUnn+/nMZzXl3j0aSXt15zZRVYwxxphIS/A6AGOMMY2TFTDGGGOiwgoYY4wxUWEFjDHGmKiwAsYYY0xUWAFjjDEmKqyAMSZEIpItIp716xeR50RERSTL51iWe+w5r+Jy4/D03pj4ZAWMiVsikuN+eIayPRdinv7pikVkn4h8IiJ/F5FzRSQxit9PTjTyjrZAhZsxtUnyOgBjgngIaBPkfBowDUgEvgoz79+7r4nue5wI/Bi4DlgpIper6jq/NFe67+mVXwF/BHZ4GENNvL43Jg5ZAWPilqo+VNM5ERHgVZwC4l/Ag2HmfU+APDOBR4BLgPkiMkRV9/qk2RrOe0Saqu4CdnkZQ028vjcmPtkjMtNQ3QtcDHwKXKURmPNIVfcAPwSyge7A3b7nA7UziOMqEVnmPmorEpFtIjJHRCa714xx0/UEetb0aM/dzxaRTu7juh0iUi4iV7vngz6mEpGBIvKGiBwUkXwR+UBEJgW47h43nzEBzh3XpuPGfpW7u9kn9pxg98Y9niAiN4nIxyKS58b1sYj8RESO+/zxuQftReQpEdnlPsb8WkSuCfR9m/hlNRjT4IjIj4BfA7uB76tqQaTyVtUKEflfYAzwIxGZWkvh9X84j64249SoDgOdgaE4NaFXgBycR3J3uGl8a2af+eXXFlgB5AGzgQpgTwih9wKWA18CM90YJgPvishlqvpKCHnU5PfABcC3gIeBQ+7xQzWmqPYicBmwDfg7oMCFwOPACODyAGnaAEuBEpzaaQrOvXxGRCpU9fk6fycmtlTVNtsazAacARS625l1SK/Or33Qa1KAUvfaXj7Hs/3TAgeA7UBagHza++3nADm1xQa8ACQFOP+cez7L51iWT7o/+10/xP0+coFWPsfvca8fE+A9KvN7rrb39jsf6N78yE3zCdDC53g6sNI9d1kN9+DvQKLP8UFAGbDa699B20Lf7BGZaTBEpBvwBpAKXK+qH0bjfVS1GKfgAOgQQpJSoDxAPvvr8PYlwM9UtSzMdIdxHhv6vv9K4GWcGsGFdYilvq51X3+pqnk+ceUDv3B3rw+QrgCYpqrlPmlW49RqThCRFlGK10SYFTCmQRCRNOBNnEc/96vqy9F+S/e1tradl3H+618tIveLyLdFpHU93jdHfToWhOETVT0a4Hi2+3pa3UOqs9NxHvFlBzi3GKdQDhTXelU9EuD4Nvc1IyLRmaizAsbEPbfH2PM4H1hvAP8T5fdLxWkLAdhXy+VT3S0P+CXwLrBfRN4Ukb51ePvddUgDNbfTVOZXn0KvrloDB1W1xP+EW0PbT+C4amrbqazVRWWckok8K2BMQ/B7nB5jXwBXqPtQPopG4HSA2aOqOcEuVNVyVX1IVb8FZAI/AF4Hvg+8JyIpYb53Xb+3zBqOd3JfD/scq3BfA3XyCTbuKFyHgbYi0sz/hIgkAe2BQDUV00hYAWPimoj8EPgNsBenx1h+lN8vgeoa0j/CSauqe1V1tqpeCiwE+gAn+VxSTvT++z5dRFoGOD7Gff3U51iu+9o9wPVDasi/sj0knPg/xfmMGRXg3Cg3r0/CyM80MFbAmLglImcAz+I0fF+oqlui/H4dgVk4H8pbgT/Ucn2KiJwd4Hgzqh+x+XahPgB0EJHmEQn4WK2B3/rFMQSnG/BhnFpVpY/c12vcmkTl9d398/BR2emhRxgxPeO+3u+2oVW+TxrOjAQAT4eRn2lgbByMiUvuf+OVPcY+BiYFGjToI0dVnwsj/3vcLxOonipmBJCM8wF8eQi9wJoDH4jIBmAVsMWNdyJwAvAfVf3G5/oFOONj3hORJUAx8Lmq/jfUuINYAlwvImfi9LaqHAeTAEzxbTRX1Q/d9x8FfCQiC3EesX0PmEPgms0C4C7gbyLyb+AocEhVH60pIFX9h4icD1wKfC0ib+A8ArwAZ9zOKzHorGG85HU/adtsC7Rx7PiOULbsEPP1T1eM09i8Cvgb8G0goYa02fiM9QCaAT/HadjfChThdApYAdwEJPulTweewBk3U4bfeJPavg+Cj4N5DqdQexPnEVgBTkFzTg15tXG/373uPfgKuJEaxsG4aaYB37jXKz5jevzvjc/xBOBmnHEvBe62Crgl0H0Odg8Cff+2xfcm7g/OGGOMiShrgzHGGBMVVsAYY4yJCitgjDHGRIUVMMYYY6LCuin7aN++vWZlZdUpbX5+Punp6ZENKAIsrvBYXOGxuMLTGONatWrVflUNPCms193Y4mkbPHiw1tWiRYvqnDaaLK7wWFzhsbjC0xjjAlZqDZ+p9ojMGGNMVFgBY4wxJiqsgDHGGBMVVsAYY4yJCutFZhqV/OIyZi7ZyEvLt5BbUEpGWjOuGN6TKaP6kJ5iv+5esZ9L02Q/WdNo5BeXceHjS9lyoIDiMmdNrYMFpcxcvIn3vtrN6zefbR9mHrCfS9NlP1XTaMxcsvGYD7FKxWUVbN6fz2/e+IqLB3cLmsfqA+Ukb6htlv7Ya8hx/WvVdjbvz6e0/NiJdYvLKthyoICZSzYybeKAaIZpPGIFjGk0Xlq+5bjCpVJpuTL70x3M/nRH7Rl9/GGEI4uQRhhXcVkFL63YagVMI2WN/KbRyC0o9ToEUwe5BSVeh2CixGowptHISGvGwSCFTFKCMDSrbY3nAQ4dyqVNm4xIh1ZvDTmuj3MOUlZR87pTGWnJkQ7LxAkrYEyjccXwnjy2cCPlARbRS0lKYMro3rU+isnOzmbMmGHRCrHOGnJc0+etZebiTQEfX6YkJXDFsB7RCs94zB6RmUbjvJM611i49GyXxpRRfTyIykwZ1Yee7dJISTr+46Zrm+b2c2nErIAxjcbf3t9c9XVSgiACbdOTmTK6t3WF9VB6ShKv33w2U0b3pm36sY/DhvdpZz+XRsx+sqZR2Lgvj9c/3V61/88bh9Xa3mJiJz0liWkTBzBt4gAWrdnLNc99DMBrK7dzy9i+dGnT3OMITTRYDcY0Cg/NX09lO/LIfu2tcIljYwZ04LQebQAoKa/gkYUbPI7IRIsVMKbBW7P7CG99sbNq/85JNqYinokIP/P5Gb22chtbDxR4GJGJFitgTIM3Y946Ktv2J5zQkVO7t/E2IFOrs/q048xeTi2zrEJ5eMF6jyMy0WAFjGnQvtpxmDlf76nanzqxv4fRmFCJyDE1zdc/3c7GfXkeRmSiIeYFjIgMEpEFIlIgIjtF5F4RSQwjfYKIrBQRFZHvBjh/voh8KSJFIrJaRCZH9jsw8WT6vHVVX593cidO7NLaw2hMOM7o1ZaR/doDUKHw8HyrxTQ2MS1gRCQDmA8ocD5wL3An8PswsrkeCDhjoYiMAP4NLALOBd4G/ikik+oRtolTq7bksnDNXgBE4I4JVntpaHxrMf/9Yidrdx/1MBoTabGuwdwENAcuUtV5qvokTuEyTURa1ZbYLaD+D/ifGi75DbBEVW9T1UWqehfwHvDbyIRv4skMn9rL+d/qQv/Mlh5GY+ri1O5tmHBCRwBUj/2ZmoYv1gXMucAcVT3ic2wWTqEzOoT09wFLgQX+J0QkBRgLvOp3ahYwXETs2UkjsmLTAT5wp4lPTBBut9pLg+Xbbvbe17v5asdhD6MxkRTrAmYgsMb3gKpuBQrcczUSkVOAa4Gf1XBJH6CZf/7ANzjfp30CNRKqyvS51f/pXnRaV3q1T/cwIlMfJ3Zpzbkndaran261mEYj1gVMBnAowPFc91wwjwCPqmpNo7Iq0/vnn+t33jRwH2zYz0c5BwFolijcNr6fxxGZ+po6sT8iztcL1+zlk625wROYBqFBTBUjIj8EBgDfi0LeNwI3AmRmZpKdnV2nfPLy8uqcNpoaW1yqyn0riqr2R3RJZOMXH7HR47iirSnEdWanRFbsKgfgN6+s4K6hdZ8+pincr0iKVlyxLmBygUBtIRlU1zSOISLNgD8DfwISRKQNUNkhIF1EWqrqUZ/0/vlX1lwC5q+qTwFPAQwZMkTHjBkT2nfix5m2vG5po6mxxbXgmz1sOrwSgOSkBP5wxSg6t47cPFaN7X5FWyTj6nFiHhNnLKG8Qvn6QAWpPU5mWO92nscVSU0trlg/IluDX1uLiHQH0ji+7aRSOk635Ok4hUQu8Ll7bhbwqfv1RqDUP393vwKwB7sNXEWF8qBP28vlZ/aIaOFivNW7QwsuOq1r1f70uevQAMsvmIYj1gXMu8A5IuLbn3QyUAgsriFNHk7vMN/tR+65u4HLAVS1GGf8yyV+6ScDy1XVuqY0cHO+3s3qXU4HxNRmCfxkjK0j0tjcNr4fSQlOY8xHOQeregqahinWBcyTQDEwW0QmuO0f9wDTfbsui8gGEXkaQFXLVDXbdwNWuJd+qaof+uR/HzBGRB4SkTEi8gBwHs6ATtOAlVcoM+ZX116uOiuLji1TPYzIREP3tmlMHtq9av9Bq8U0aDEtYFQ1FxgPJAL/xRlkOQP4nd+lSe414eb/AXAxMAGYA3wfuExV59YjbBMH3vpiJ+v2OHNVpScn2iqIjdit4/qS7K5++dm2Q1WzNZiGJ+a9yFR1NTCulmuyajmfA0gN594A3qhjeCYOlZVX8JDPPFXXjeh13MqIpvHo3Lo5l5/Zg2eX5gDOuJixAzqSkBDwT97EMZtN2cS92Z/uYPP+fABapSZx3cjeHkdkou0nY/qQ2sz5ePp65xHmfL3b44hMXVgBY+JaSVkFf/VZK+SGkb1p3byZhxGZWOjYMpWrhmdV7c+Yv47yCmuLaWisgDFx7bVV29ieWwhARlozrhnRy+OITKxMGd2H9GSnKXbdnrxjVi01DYMVMCZuFZWW88iC6pmBbhrdhxYpDWLyCRMBbdOTudbnH4qH5q+nrLzCw4hMuKyAMXHrnx9tZfcRZ1qY9i1SuNLnkYlpGq4f2ZtWqc4/FZv35/P6pzs8jsiEwwoYE5cKS8p5bFH1DGO3jO1D8+Swe66bBq5182bc4NOp468L11NSZrWYhsIKGBOXXliew/68YgA6t07lR2f08DYg45lrRvQiI83p2LHtYCGvrdrmcUQmVFbAmLiTV1zGk4uray+3jutLajOrvTRVLVKSuGl09cDaRxduoKi03MOITKisgDFx59kPNpNbUApAt4zmXDK4ey0pTGN35fAs2rdIAWDX4SL++dFWjyMyobACxsSVw4Wl/O39TVX7t43vVzVtiGm6micncrPP5KaPLdpIYYnVYuKd/eWauPL0+5s4UlQGQK/26cdM326atsvO7EGnVs4Ep/vzinlxRY63AZlaWQFj4sbB/BKe/mBz1f4dE/qRlGi/osaR2iyRW8f1rdp/InsjecVlHkZkamN/vSZuzFyykXz3sUf/zBZ895QuHkdk4s2lQ7rTLcNZZC63oJTnlm6uJYXxkhUwJi7sPVrE88tyqvanTuhPos2ea/wkJyVw2/h+VftPLdnE4cJSDyMywVgBY+LCE9kbKSp1BtAN6tyKc07s5HFEJl5ddFpXerVPB+BIURlP+3QKMfEl5gWMiAwSkQUiUiAiO0XkXhEJOshBRE4Ukffc64tFZKuI/F1EOvtd95yIaIBtYHS/K1Mfuw4X8vKH1d1O75zU39b+MDVKSkzgjgnVtZhnluZwML/Ew4hMTWJawIhIBjAfUOB8nKWM78RZ2TKY1sBm4GfAOTgrYE4A3hER/9kP1wDD/bacyHwHJhoeXbihavqPU7u3YdzAjh5HZOLdd0/pQr+OLQBnYO7MJRtrSWG8EOupaW8CmgMXqeoRYJ6ItALuEZEH3GPHUdVlwDKfQ9kish2YC5wCfOJzLl9VV0QnfBNp2w4W8OrK6qk/7pzUHxGrvZjgEhOEaRP785OXnT/955flcN2IXnRsmepxZMZXrB+RnQvM8StIZuEUOqPDzOuA+2pr5zZgf12wntJyZyGpM7LaMqJve48jMg3FOSd2YlDnVgAUlVbwRLbVYuJNrAuYgTiPsKqo6lagwD0XlIgkiEiyiAwA/gh8DHzkd9kgETnittV8ICLhFlwmRjbvz2e2z/Tr06z2YsKQ4NZiKr384VZ2HS70MCLjT1RjtwypiJQCd6nqQ37HtwMvqOrdtaR/D6cNBmAVcJ6q7vU5fztQAqwGOuC07wwGRqiqf0FUmeZG4EaAzMzMwbNmzarLt0ZeXh4tWrSoU9poiue4Xt6YxPJdzriXE9slcNfQ5h5HFd/3y+I6nqpy34oiNh122vDGdU/iyhNTPI+rJo0xrrFjx65S1SEBT6pqzDagFLgjwPHtwB9CSN8POBO4AqcmtApIDXJ9Gk7ngDdCiW/w4MFaV4sWLapz2miK17he+s8CzfrlW9rzF862astBr0NS1fi9XxZXzRav3Vv1e9T37rd164H8uIgrkMYYF7BSa/hMjfUjslycHmH+MtxzQanqelX9UFVfwqnJnAZcFuT6AuAd4PS6hWui5fUNJVRWnscN7MjpPTK8Dcg0WCP7teeMrLYAlJYrjyxc73FEplJYBYzb/tFVRPq4XY7DtQa/thYR6Y5T01gTMEUNVHULcBDoXdul7mbixNc7D7NyT/VMuL7P0Y0Jl4gwbVL179C/P9nB7nxb9TIe1FrAuIMc/yQiq4A8YCuwDtgvIntF5A0RuUJEQnmA/i5wjoi09Dk2GSgEFocTuNvQ3w7nEVhN1zQHvoPzKM3EiRnz1lV9fc6JmZzUNVCl1pjQDevdjrP7tgOgvEJ5c4MNvIwHNRYwInK2iCwCvsDpQrwYuAH4Ps7jqcnAA0ARMAPYKSK/EZFgLUVPAsXAbBGZ4Daw3wNMV5+uyyKyQUSe9tn/i4j8UUQuFJGxInIzMAfYiNPNGRFpLSLvi8gUERkvIpOBRUAX4A9h3hcTJZ9tO8T8b5x+GSIw1WovJkKmTRxQ9fWKXeWs33PUw2gMBB9oORv4K/BjVd0eLBN3qpcJwB3uofsCXaequSIyHngU+C9wCKdwuidAXL7Tx6wEforT2ysVpxb1b+B+Vc13rykG9gG/BjriFHzLgdGqujJY/CZ2Hpy7turr757ShYGdWnkYjWlMBvfMYOyADixauw8FZsxfx+OXD/Y6rCYtWAHTU1WLQslEVctxahRzRCToUFpVXQ2Mq+WaLL/9Wbg1lSBpioCLQonXeOPjnIO8v34/AALHzCdlTCRMmziARWv3AfDOl7v5eudhTuxij2C9UuMjslALl0ilM42bqvKXOdW1l7O6JNGnQ/yNBzAN28ndWnPOiZlV+zPmWY8yL4XdTVlEkkTkZhF5TUT+JSK3BJhw0phjLNt4gA83HwQgKUE4v28zjyMyjdXUif2pnA9i/jd7+GzbIU/jacrqMg7mYeBKIBtnmpafA49FMCbTyKjqMW0vlwzpTsc0W4rIRMfATq04o1N1E+50n16LJraC9SKraXDixcAkVX1MVR8AbgYuiUZwpnHIXruPT7Y6/0UmJybwU5911Y2Jhgv6JlO5pNCSdfv4OOegtwE1UcH+jXxXRP4mIh38ju8GJgKIMzPhWGBXlOIzDZyqHvMf5I/O6E6XNt7POWYat84tErjgtK5V+741aBM7wQqYATizHH8jIneJSOVD81uBR0RkL870LpcBP4lumKahmrt6D1/uOAxASlICt4y12ouJjdvH9yPJrcas2HSQZRv2exxR0xOsF9khVb0dGIlTS1ktIt9X1feBLJxxL6OAXqq6JBbBmoalokKZPre69nLl8J50bGULQpnY6NkunUuGdKva/8vctZWT4JoYqbWlVVW/UdXzcAZR/klE5gF9VfULdyuOepSmQXr7y12sdUdTpyUnctPoPh5HZJqaW8f1IznR+Zj7ZOshstft8ziipiWUuciai0hrVX0bOBlnPrHFIvKoiLSNeoSmQSorr2DG/OrayzVnZ9GuRYqHEZmmqGub5vzojO5V+9PnrrNaTAwF60XWV0SWAvnAQRFZjzPtynRgEM4sAGtE5KfuVDHGVHnzs51s2ufM4tMyJYkbRtY26bUx0XHL2L6kJDkfdV/uOMzc1Xs8jqjpCFaDeQHIAToBbYDncSapTFXVfap6E05vsguBL6MdqGk4SssreHhB9Qjq60f2pk1asocRmaasY6tUrhzes2p/xrx1VFRYLSYWghUwg4DnVHWvqh4FHgdaAlWtZqr6uaqOA/4numGahuRfq7az9WABAG3SmnHtiCxvAzJN3k2j+5CW7DxoWbP7KG9/aSMrYiFYATMXp1H/ByJyHvAMzvT4G/0vVNXXoxSfaWCKy8p5xKf2MmVUH1qm2rQwxlvtWqRwzdlZVfsz5q+jrNwWJYu2YAXMtTiFzC+B+3FWj5yg1kJmgpj10TZ2HnbmO23fIpmrzupZSwpjYuOGkb1pmeJMm7hpXz5vfrbT44gav2DjYPJU9ZeqOlRVv6Wq17rLFBsTUFFpOY8t2lC17zyWsHlQTXxok5bMdSN7Ve0/vGA9pVaLiapgvcikpnPB1JZORAaJyAIRKRCRnSJyb2290Nxlm99zry8Wka0i8ncR6Rzg2vNF5EsRKRKR1e7KliYGXlqxhb1HnWFRma1SuGKY1V5MfLl2RC/apDmPbLceLODfq4KupWjqKdgjsnUicr2IpIeSkYgMFpEXcB6p1XRNBjAfUOB84F7gTuD3tWTfGtgM/Axnuebf4cwk8I7vUgEiMgJnpctFwLnA28A/RWRSKN+Dqbv84jIez65unrt1bF9Sm1nvdRNfWqU248ZR1V3m/7pgPcVl5R5G1LgFe37xC5wP/odFZC6wDPgK2I+zPHEboBcwGPg20B34O05ngJrcBDQHLlLVI8A8EWkF3CMiD7jHjqOqy9z3r5QtIttx2ohOAT5xj/8GWKKqt7n7i0TkROC37rUmSp5blsPB/BLAGdx26dDutaQwxhtXn5XFMx9sZn9eCTsPF/HKx9u4cniW12E1SsHaYGbjfHhfAOThTBXzNrAC+BRnPZiZOAXMTJwlln+qqsFGMZ0LzPErSGbhFDqjw4z9gPuaDCAiKThzpr3qd90sYLiI2LqpUXKkqJSnlmyq2r9tfF9Skqz2YuJTWnLSMdMWPbpwA0WlVouJhqBTxahjnqr+WFW7Al2BIcAI4ASgtaqeraoPqureEN5vILDG7z224szaPLC2xCKSICLJIjIA+CPwMc6iZwB9gGb++QPf4Hyf/UOIz9TB0+9v5nBhKQA926Vx0endaklhjLeuGNaTzFbO1EV7jxbz0grrvxQNEstexyJSCtylqg/5Hd8OvKCqd9eS/j2cNhiAVcB5lQWbiJwNfACcpqqf+aTpC6wHzlHV4x6TiciNwI0AmZmZg2fNmlWn7y0vL48WLeJvjflox5VXoty1pIDCMmf/xlNSOKtL7T3Hmur9qiuLKzyhxLVgaykvrnYe67ZMhj+PSiM1qU59myIalxfqE9fYsWNXqeqQQOcaWh/SnwJtgX7Ar3EWRTtbVYvqmqGqPgU8BTBkyBAdM2ZMnfLJzs6mrmmjKdpx/em9NRSWOY37fTu24Bc/HEViQu1/pE31ftWVxRWeUOIaXlbOwr8sZsehQo6WwMbE7twyJrrrFTXk+1UXsV4YPRenR5i/DPdcUKq6XlU/VNWXcGoyp+EseFaZNwHyz/A7byJkf14xzy3Nqdq/Y0K/kAoXY+JBSlLiMct3P7VkE0eKSj2MqPGJdQGzBr+2FhHpDqRxfNtJUO6gz4NAZZ/DjUCpf/7ufgWwDhNRT2ZvpNBtHB3YqSXnnXTcsCRj4toPBnejZ7s0AA4XlvLMB5s9jqhxiXUB8y5wjoi09Dk2GSgEFoeTkdvQ3w5nfAzuwmeLgEv8Lp0MLFfVw3UN2hxvz5EiXvRpGJ02sT8JVnsxDUyzxARuH9+vav/p9zdzqKDEw4gal5AKGBE5OULv9yTOGJrZIjLBbWC/B5ju23VZRDaIyNM++38RkT+KyIUiMlZEbgbm4NRafFvl7wPGiMhDIjJGRB4AzsMZ0Gki6LFFGyguc6bZOKVbayYOyvQ4ImPq5vxTu9KngzOe/Ghx2TFd7k39hFqD+VxEPhaRn4hIm7q+marmAuOBROC/OAM5Z+CMzPeV5F5TaSUwEngaZyzObTgj9oepar5P/h8AF+OM8p8DfB+4LFDvMVN323ML+OdHW6v2p03sTx1nFjLGc4kJwtSJ1aMYnl2aw/48Wwk+EkLtRTYOuAZ4AHhQRN7EGbE/P9zZlVV1tZtfsGuy/PZncWxNJVjaN4A3wonJhOfRhRsoLXd+7IN7ZjC6fwePIzKmfs47qTMDO21gze6jFJaW82T2Rn793UFeh9XghVSDUdVsVb0KZ3XLW3EGXM4BtojIfSLSJ2gGptHI2Z/Paz4TBN45yWovpuFLSBCm+dRiXlyxhT1H6jz6wbjCauRX1XxVfUZVRwEDcJZUvhtnYszFInJhFGI0ceSvC9ZT7i43O7x3O87q097jiIyJjImDMjmlmzPKobis4pilJ0zdhN2LTESyROQenBrMcOAdnJHwe4BXRGRGRCM0cWPD3qO88dmOqv07J9nsO6bxEDm2FvPPj7ayPbfAw4gavlB7kaWJyJUisgjYAFwO/A3ooarfU9WnVfVSYApwXfTCNV6aMX89buWF0f07MCSrrbcBGRNho/t3YHBPZ2x2abny6EKrxdRHqDWYPcATwHacZZP7qer9qrrL77qPqZ7l2DQi3+w6wttfVP+4ff/TM6axEBHu9Pndfm3VdrYcyA+SwgQTagHzc6CLO6tydk0XqepXqtqrpvOm4Zoxr3oihImDMvlW9zr3Vjcmrp3Vtz3De7cDoLxCeXjBeo8jarhCLWA6AAFXthSRziLy28iFZOLNF9sPMXd19TI/VnsxjZ1v++Ibn+5gw948D6NpuCoEK6cAACAASURBVEItYH4H1LTIRxeOHyhpGpHpPrWX75zSmRM6t/IwGmOib0hW26rxXRUKD823qQzrItQCRoCaBlR2w2YqbrRWbTlI9tp9ACQITJ3Qr5YUxjQOvjX1t77YxTe7Aq7oboKocSS/iFwFXOXuKvCEiPjf4VTgZGy9+0brwbnV/7mdf2pX+nZsGeRqYxqPb3Vvw4QTMpn/jfN4eMa8dTx1ZcB1tUwNgtVgCnB6hB3AqcEc9tmv3DbjTB9zY3TDNF5YtnE/yzY6nQITE+SYWWeNaQp8azFzV+/hy+02KXs4aqzBqOprwGsAIvIscJ+q2jSjTYSqMt2n9nLx6d3Iah+wn4cxjdagLq34zsmdeftLp4v+g/PW8tw1Z3gcVcMR6lxk11jh0rQsWb+flVucprVmicJPx0d3KVlj4tXUif2oXOooe+0+Vm056G1ADUiwNpgHgL+q6nb362BUVX8R2dCMV1SVB+eurdr/4dAedMtI8zAiY7zTt2NLzj+1K69/6kyT9ODcdfzjhmEeR9UwBJuu/xLgZZzR+/6rRPpTwAqYRmL+N3v5wn3WnJyUwC1jrfZimrbbx/fjP5/vpLxCWbbxAMs3HmB4n3ZehxX3anxEpqq9VPVzn6+Dbb1DfUMRGSQiC0SkQER2isi9IpJYS5qhIvKsu9JlgYisFZHfiUiq33X3iIgG2L4danxNXUWFHjPu5cfDetKpdWqQFMY0flnt07n49OqhgNPnrSXMpbCapLBnU64PEckA5uPUeM7HWcr4TpyVLYOZDPQB/oSzBPJjwDScGpa/wzizPPtuyyMQfpPw7le7q/r7N2+WyE/G2FI/xgD8dHxfmiU6jTEf5+SyZP1+jyOKf8HaYM4LJyNVfSeEy24CmgMXqeoRYJ6ItALuEZEH3GOB/FFVfX+a2SJSBMwUkZ6qusXnXJmqrggnduMor1Bm+IxYvuqsLNq3SPEwImPiR7eMNCYP7c5LK5zlwqfPXcuofu1twb0ggrXBvIVT0wjl7ikQ9DGX61xgjl9BMgunZjIa+G/AzI8tXCp96r52AbYEOG/C9J/Pq+dcapGSxJRRIT/5NKZJuHVsP15duZ2Ssgo+336Y+d/sZeKgTK/DilvBHpH1Anq7r7VtoX4SDQTW+B5Q1a04gzoHhhM4zqOvCmCj3/E2IrJfREpF5FMRuSjMfJuksvIKHp5fPWvstSN6kZGe7GFExsSfTq1TueLMnlX70+eto6LC2mJqIrFsqBKRUuAuVX3I7/h24AVVvTvEfDoBXwDvqOrVPsevADri1G5a4iyAdh7wA1WdXUNeN+LORJCZmTl41qxZ4X5bAOTl5dGiRYs6pY2mUONasr2UZ74qASC9GTwwKo30ZtGr+jf0+xVrFld4ohnX4WLlriUFlJQ7+7ecmsLQTsEeBsUmrvqoT1xjx45dpaqB59BR1YAbkOb7dW1bTfn45VkK3BHg+HbgDyHmkQwsATYBGbVcKzgN/J+FkvfgwYO1rhYtWlTntNEUSlzFpeV61v0LtOcv3tKev3hLH124Pi7i8oLFFZ6mGtf973xT9fcy/sFsLSuviIu46qo+cQErtYbP1GCPyI6KSOWcCHnA0Vq2UOQCrQMczyCEGZnFaU17ATgROE9Vg6Zxv/nZwCm1dYVuyl5ZuY0dhwoBaJuezNVnZXkbkDFxbsqo3rRIcWotG/bm8d/Pd3ocUXwKVq+7lur2jWupebr+cKzBr61FRLrj1ILWBExxrIdwujdPVNVQrgcnbntIWoOi0nIeXVjd9vKT0X1ITwmtum9MU5WRnsy1I3rxV3e1y4fmr+O7p3QmKTGmIz/iXrDJLp/3+fq5CL3fu8BdItJSVStrPZOBQmBxsIQi8ivgVuBSVf0glDdzazw/AD5X1fK6h914vfzhVvYcKQagQ8sUrhjWs5YUxhiA60b04rmlmzlSVEbOgQJmf7KDS4d29zqsuBJWcSsibURkhIhc4r6GuzD7k0AxMFtEJrgN7PcA09Wn67I7Yv9pn/3LgD/gPB7bISLDfLYOPtctFpHbRGSSiFwIvA2c6b6H8VNQUsYT2Ruq9m8d25fmyfYk0ZhQtG7ejCmjqwciP7xgPSVlFR5GFH9CKmBEJElE/oTTGL8EeMV93S4iD4hIs1DycdtMxuOMmfkvzgj+GRy/5HISx46rmeS+Xo3TaO+7fcfnug3AHcCbwEs4Pcm+o6r/CSW+pub5ZVvYn+f0HOvSOpUfnmH/fRkTjqvPyqKt251/x6FCXlm5zeOI4kuoD9un43TlvRen0XwvTnfgHwC/xlnZ8rZQMlLV1cC4Wq7J8tu/GqdwqS3v60KJwcDRolJmLqkeQnTruH6kJFntxZhwpKckcdPo3vzhHadJ+LGFG7hkcDdSm9nfEoT+iOzHwN2q+gdVXaOqB93X/8MpYH4cvRBNNDy7NIdDBaUA9GibxiVDutWSwhgTyI+HZdGhpTOl0u4jRfzjw60eRxQ/Qi1gKoCvazj3FdZLq0E5XFDK396vXj/utvH9aGa9X4ypk+bJidziMyns49kbKCgp8zCi+BHqp8qLwPU1nLsBp73DNBB/e38TR4ucP4DeHdK54NQuHkdkTMP2ozN70MVd1mJ/XgkvLLfpESFIASMiN1duQA4wXES+FpH7RWSq+7oap5eW/3xgJk4dyCvmmaWbq/bvmNDf+u4bU08pSYncOq5f1f6TizdytKjUw4jiQ7BG/kcDHOsCnBDg+HTg4YhEZKJq5pJNFLiTKA3IbMl3T+7scUTGNA6XDOnGE4s3sO1gIYcKSnl2aQ63je9Xe8JGLNiKlglhbNZlogHYe6SIF5bnVO1PndifhARby8KYSGiWmMDt4/tX7f/t/U0cLmjatRh7NtKEPJ69kaJSZyDYSV1bcc6Jto6FMZF0wald6N0+HYCjRWXHdKZpisIdyd9NRMaJyHn+W7QCNJGx81DhMd0n75w4wFbiMybCkhITuGNidS3m2aWbOZBX7GFE3gp1JH9LEXkXZ+XIeTirXb6FMxq/cjNx7JGFGygpd2ovp/Vow5gBHWpJYYypi++e3JkBmS0ByC8pZ+aSpluLCbUGcz/QAxiJs8bKhcAY4GlgMzAsGsGZyNh6oIDXfKawsNqLMdGTkCBMnVjduP/C8hz2Hi3yLiAPhVrAnAf8H/Chu79TVZeo6o04837dFY3gTGT8deF6ytxlXc/s1Zaz+7bzOCJjGrdzTuzEiV1aAVBUWsHji5rmSI5QC5hMYJs75X0+0Nbn3DtUT0Zp4syuvApmf7K9av/OSVZ7MSbaRIQ7J1W3xfzjw63sdBf1a0pCLWC2Ae3dr9cD3/U5dybQNOt/DcCbG0twKy+M7NeeM3q1DZ7AGBMRYwd05NTuzoomJeUVPLpoQy0pGp9QC5h5wAT36xnALSKyTEQWAffhrNNi4sza3Uf5cFf1OmvTfHq3GGOiS0T42aQBVfuvfryNvQVNa72YUAuYX+Cu2aKqL+JM078ZyMVZZfKXUYnO1MuMeeuqZiEdP7Ajp/XI8DQeY5qas/u2q3pqUFah/Gdj0xp4GVIBo6oFqrrfZ/91Vb1cVS9S1SdUNeRiWUQGicgCESkQkZ0icq+IBJ0JQESGisiz7kqXBSKyVkR+JyKpAa49W0Q+FJEiEdksIiGtU9PYfLXjMO99vbtqf6rVXoyJORHhTp+/vaU7yti0L8/DiGIr3IGWA0TkChG5S0QuF5EBtac6Jn0GMB9nev/zcRYwuxNnZctgJgN9gD/h9Gh7DJgGvOyXf19gDk7t6jxgJjBdRGqaCbrRmj5vXdXX557UiZO6tvYwGmOarjN7t2NkP6cJW4GH5q/3NqAYCmlFSxFpBfwN59FYApAHtAAqRGQ2cL2qHgkhq5uA5sBF7vXz3LzvEZEHguTxR98aFJAtIkXATBHpqaqVc2PfBewErlDVMmChiPQAficiT6tqk1i35pOtuSxcsxdwBi1Z7cUYb02b2J/31zsfYf/9Yie3jO3LgE4tPY4q+kKtwTyO0xX5SiBdVVsB6cBVwET3fCjOBeb4FSSzcAqd0TUl8itcKn3qvvouZnIuMNstXHzz7wacFGKMDd4Mn9rLmZ0T6Z/Z+H+RjYlnp/XIYPzAjgCowkPz19WSonEItYA5H7hLVf+hqoUAqlqoqi8DP3fPh2IgsMb3gKpuBQrcc+EYjrPS5kYAEUkHuvvnD3zj896N3oebDlT9p5QgcEHfZI8jMsbAsU8S3v1qN1/tOOxhNLER0iMynEdiu2o4txNn8GUoMoBDAY7nuudCIiKdgF8DL6rqXvdwG/fVP/9cn/cOlNeNwI0AmZmZZGdnhxrGMfLy8uqcNlJUlfs/qh6SdHaXJFpogedxBRIP9ysQiys8Fld4Tm2nfHbAGej8P7OWM3Xwcf2UPBGt+xVqAfMY8DMRWVhZgwEQkTTgZ4T+iKzeRCQZeBWn0Jta3/xU9SngKYAhQ4bomDFj6pRPdnY2dU0bKe+v38e6OR8BkJQg/OHykWz84iPP4wokHu5XIBZXeCyu8Ow4upDPlxWiCp/vK6dV729xehwMH4jW/Qq2ZPIDlRvQCugHbBORf4rIwyLyT2Ar0BcI9SF/LhCoO1MG1TWNGokzx8kLwInAearqm6ay5uKff+VPr9b8GzJV5cG51c91Lx3ane5t0zyMyBjjr2vLBL53SnWzsW97aWMUrAZzid9+qbv5zpx81H39AaFNeLkGv7YQEekOpHF820kgD+G090xUVf+2nHwR2eafv89+KPk3WIvW7uWzbU4Zm5yUwE/H9fU4ImNMIHdM6MdbX+ykQuH99fv5cNMBzuzdOCegDbZkcq8wtt4hvt+7wDki4lvjmQwUAouDJRSRX+HMGnCFqn4QJP8L/QZuTsaZS+2rEGNscPxrL5ed0YPOrZt7GJExpia9O7TgotO7Ve0/OG8djXUERayXTH4SKAZmi8gEt4H9HmC6b9dld8T+0z77lwF/wHk8tkNEhvlsvitn/RmnS/KLIjJWRH4OTAHubcxjYOZ8vZuvdzq3L7VZAjeP7eNxRMaYYG4f34+kBKex/6PNB1m64YDHEUVHyAWMiPQWkSdE5EsR2eG+Pi4iodZecNtMxgOJOKtg/h5n8szf+V2a5F5TqXI5gKuB5X7bd3zy3wB8G6dd6F3gZuBOVf17qDE2NOUVesyo/auGZ9GxZXz0TDHGBNa9bRqXDu1etf+XuWsbZS0m1JH8g4FFONPyvwXswVkj5gfA5SIyVlU/CSUvVV0NjKvlmiy//atxCpdQ8v8AOCOUaxuDt77Yybo9ztxG6cmJTBlttRdjGoJbx/blXyu3U1JewWfbDrFo7V7GDcz0OqyICrUG8xeckfNZqnqtqv5KVa8FernH/xKtAE3NysoreNhnXqNrR/SibboNrDSmIejSpjmXndmjav/BuY2vLSbUAuYM4AFVLfA96O7/BWfRMRNjr3+6g037nTGuLVOTuH5EyE8rjTFx4OaxfUht5nwMf73zCHN8ZkBvDEItYAqBmvrRtcVWtIy5krIK/rqwuvZyw8jetE5r5mFExphwdWyZypXDs6r2p89bR3lF46nFhFrAvA38UURG+B509+/HabA3MfTaqm1sO+hMqpCR1oxrzs7yNiBjTJ1MGdWb9GSnT9O6PXm89cVOjyOKnFALmGnAJmCxiOwSkc9FZBfO2JXNOGu6mBgpKi3n0YXV63tPGd2HlqlWezGmIWrXIoVrzu5Vtf/w/PWUlTeOpZVDXdHygKqOwOkS/Diw1H09V1VHqGrj7MQdp2Z9tJVdh52nku1bpHDl8J4eR2SMqY8bRvamZarTqXfT/nze+Kxx1GJqLWBEJFVE1onIt1X1PVW9T1Vvdl/nxiJIU62wpJxHF22s2r95TB/SkkOds9QYE49apzXjhpHVnXQeXrCO0kZQi6m1gFHVIpyp8Bv+d9sIvLgih/15xQB0apV6TDdHY0zDdc3ZWbRxO+psO1jIayu3exxR/YXaBvMycE00AzG1yysu44ns6trLreP6ktosMUgKY0xD0TK1GTf5DJR+ZOF6ikrLPYyo/kJ9trIVuFREPsaZgmUP4NuXTlX1iUgHZ4713NLN5BaUAtAtozmXDuleSwpjTENy5fCe/P39TezPK2HX4SJmfbSVq306ADQ0oRYwD7qvnYHBAc4rYAVMFB0uLOWpJZuq9m8b34/kpFjPVWqMiaa05CRuHtOXe99aDcBj2RuZPLQHzZMb5pOKUHuRJdSyNczvvgF5+v1NHCkqA6BX+3QuOq2rxxEZY6LhsjN70KmVM2HtvqPFvLgix9uA6sH+BW4AcvNLeGZpTtX+7eP7kZRoPzpjGqPUZonc4rNg4JOLN5FXXOZhRHUXznT9ySJyo4j8XUTedl9vEBGbXTHKZi6p/gXr17EF3/tWl1pSGGMasslDutO1jbNo4MH8Ep5fluNtQHUUUgEjIicA64HHgJOAcvf1MWCDiAyKWoRN3L6jxcf8ck2d2J9Ed6EiY0zjlJyUwO3j+1Xtz1y8kcOFpR5GVDeh1mCeAg4DfVR1mKp+X1WH4SzsdQhnpcqQiMggEVkgIgUislNE7vVb4jhQmmQR+bOIvC8ihSIScDY4EXlORDTANjDU+OLNE9kbKXS7Kg7q3Ipvn9jJ44iMMbFw0eldyWqXBsCRojKe/mCzxxGFL9QCZgjwW1Xd6nvQ3f8dMDSUTEQkA5iP0+vsfOBenHnMfl9L0jTgeqAAWFbLtWuA4X5bTijxxZtdhwt56cMtVfvTJvYnwWovxjQJSYkJ3DGhf9X+Mx9sJje/xMOIwhdqAZMD1LQObyrOOJlQ3AQ0By5S1Xmq+iRO4TJNRFrVlEhVDwFtVfUc4PVa3iNfVVf4bQ1yOYHHFm2gpMyZQOFb3dsw/oSOHkdkjIml732rC/06tgCcgdYzfYYqNAShFjC/BP5XRI5ZWExEhgH3Ab8IMZ9zgTmqesTn2CycQmd0sITa2JZ6q8W2gwW88vG2qv07J/ZHxGovxjQliQnC1InVtZjnl+Ww72ixhxGFJ9QC5tdAK2CZ33T9S93jd4vIR5VbkHwG4jzCquI+Zitwz0XCIBE5IiLFIvKBiAQtuOLVIwvXU1rulKlDszIY2a+9xxEZY7zw7RM7cUJn5wFPYWn5MdNFxTsJpWIgIs+Gk6mqBpy3TERKgbtU9SG/49uBF1T17hBiuRV4RFWP+3deRG4HSoDVQAec9p3BwAhVDVjwiciNwI0AmZmZg2fNmlVbCAHl5eXRokWLOqX1tzu/grs/KKRyYbtfnpHKwLZ1G8saybgiyeIKj8UVnsYW16d7y3j4E6fmkpQAfx7VnIzUyI2Fq8/9Gjt27CpVHRLoXEhTxdRUYMQbVX3Yd19E3gG+Bu4GLqghzVM4veQYMmSIjhkzpk7vnZ2dTV3T+rtj1qdUqLNa5dl923HTRcPqnFck44okiys8Fld4Gltco1XJ3reMz7cdoqwCVhV14H+/fbLncdUm1sPBc4HWAY5nuOciSlULgHeA0yOdd7Ss33OUNz+vXmxo2sQBHkZjjIkHIsI0n7aYVz7exvbcAg8jCk2sC5g1+LW1iEh3nG7IawKmqD/l2Jmf49pD89dT+dRy7IAODO6Z4W1Axpi4MKpfe4ZmOZ8HpeXKIws21JLCe7EuYN4FzhGRlj7HJgOFwOJIv5mINMdZ5nlVpPOOhq93HubtL3dV7VvtxRhTyanFVH8m/OuT7eTsz/cwotrFuoB5EigGZovIBLeB/R5gum/XZRHZICJP+yYUkXNF5GLgVHf/Ynfr6e63dkf6TxGR8SIyGVgEdAH+EJPvrp5mzFtf9fWkQZmc3C3Q00RjTFM1vE87zurTDoDyCuXhBetrSeGtmBYwqpoLjAcSgf/iDLKcgTMbgK8k9xpfTwCvAde5+6+521h3vxjYh9Ol+h2chvtDwGhVXRnRbyQKPtt2iPnf7AFABKZN6l9LCmNMU3Snz2fDG5/tYP2eox5GE1yoC45FjKquBsbVck1WKMf8zhcBF9UnNi9Nn7eu6uvvnNyZgZ1qnNjAGNOEDe7ZljEDOpC9dh+qTrvtY5fHZz8mW1QkDnycc5Al6/YBkCAcM/+QMcb4u9OnLebtL3exeueRIFd7xwqYOPDg3LVVX19wWlf6doy/AWLGmPhxcrfWTBqUWbXv+wQknlgB47FlG/azYtNBwJl3yHcNCGOMqYnvHGXzv9nD59sOeRhNYFbAeEhV+YtP7eXSId3o2S7dw4iMMQ3FCZ1b8d1TOlftPxiHtRgrYDyUvW4fn2x1/utITkzg1nFWezHGhO6OCf2pXCJqybp9rMw56G1AfqyA8YiqMn1u9X8cPzqjeg1uY4wJRd+OLbjgtK5V+w/Oja9ajBUwHpm7eg9f7jgMQEpSAreM7etxRMaYhuj28f1IdKsxyzcdYNmG/R5HVM0KGA9UVCgzfJ6X/nhYTzq2qmnBUGOMqVnPdulcMrhb1f6D89YRL+szWgHjgXe+2sWa3c7o27TkRG4a08fjiIwxDdlPx/cjOdH5OF+1JZfF7rg6r1kBE2PlfrWXq8/Kon2LFA8jMsY0dF3bNOeHZ3Sv2p8eJ7UYK2Bi7M3PdrBxnzMDasuUJG4c1dvjiIwxjcEtY/uSkuR8pH+x/TDzVu/xOCIrYGKqtLyCh+ZXz3563chetElL9jAiY0xjkdkqlR8P61m1P33eOioqvK3FWAETQ/9etZ2tB51V6Fo3b8a1I3p5HJExpjG5aUwf0pKdiejX7D7KO1/tqiVFdFkBEyPFZeU8srB6Bbopo3vTKrWZhxEZYxqb9i1SuPqsrKr9GfPWUe5hLSbmBYyIDBKRBSJSICI7ReReEfFf+8U/TbKI/NldUKxQRGq8YyJyvoh8KSJFIrLaXXjMc698vI0dhwoBaJeezFXDs7wNyBjTKN04qjctU5yVWDbuy+fNz3Z4FktMCxgRyQDmAwqcD9wL3Imz8FgwacD1QAGwLEj+I4B/46xkeS7wNvBPEZlU7+Droai0nEd9ai8/GdOH9JSYL8VjjGkC2qQlc93I6sfvDy9YT2l5hSexxLoGcxPQHLhIVeep6pM4hcs0EalxhS1VPQS0VdVzgNeD5P8bYImq3qaqi1T1LuA94LeR+xbC99KKLew9WgxAZqsUrvBpiDPGmEi7dkQvWjd3HsFvOVDAv1dt9ySOWBcw5wJzVNV3dZxZOIXO6GAJtZZO3SKSgrN88qt+p2YBw0XEkwXu84vLeCJ7Y9X+LWP7ktos6BNBY4ypl1apzY4ZAvHIwg0Ul5XHPI5YFzADgTW+B1R1K86jr4H1zLsP0Mw/f+AbnO/Tk2Uin1+ew4H8EsAZDDV5aPfgCYwxJgKuPiuLdunOMIgdhwp59eNtMY8h1gVMBhBoVZxc91x98yZA/rl+52PmSFEpMxdvqtr/6bi+pCRZ7cUYE33pKUn8xGcaqkcWbqCoNLa1mCbf0iwiNwI3AmRmZpKdnV2nfPLy8o5L+8aGEg4XlgLQobnQPm8j2dmbAqSOnkBxxQOLKzwWV3gsLkePcqVNinCoWNl7tJjfv7yQc7KOHx4RrbhiXcDkAoHaQjKormnUJ28C5J/hd/4YqvoU8BTAkCFDdMyYMXV68+zsbHzTHioo4aeLFlXt/+p7pzDh9G4BUkaXf1zxwuIKj8UVHour2u60HH775tcAzN0Gv/7RiON6sUYrrlg/IluDX1uLiHTH6Ybs33YSro1AqX/+7n4FENOVeJ5asomjxWUA9OmQzvmndq0lhTHGRN7kod3p0tpZDuRAfgnPL8+J2XvHuoB5FzhHRFr6HJsMFAKL65OxqhbjjH+5xO/UZGC5qh6uT/7h2J9XzHPLcqr2p07sX7UgkDHGxFJKUiK3ja9ejn3m4k0cKSqNyXvHuoB5EigGZovIBLf94x5gum/XZRHZICJP+yYUkXNF5GLgVHf/YnfzHVRyHzBGRB4SkTEi8gBwHs6Azph5MnsjBSVOY9rATi0576TOsXx7Y4w5xg8Gd6NH2zQADheW8swHm2PyvjEtYFQ1FxgPJAL/xRlkOQP4nd+lSe41vp4AXgOuc/dfc7exPvl/AFwMTADmAN8HLlPVuRH9RoLYc6SIF1dsqdqfOrE/CVZ7McZ4qFliArf71GKefn8zhwpKov6+Me9FpqqrgXG1XJMVyrEa0r4BvFGX2CLhsUUbKC5zpmU4uWtrJg3K9CoUY4ypcsFpXXk8ewMb9+VztLiMp5Zs4uffru/ww+BsNuUI2nGokFkfVQ9mmjapPyJWezHGeC8xQbhjQvV48+eW5XAgrziq72kFTAQ9unA9Je6kcoN7ZjCmfwePIzLGmGrfObkzAzs5fawKSsp5cvHGWlLUjxUwEbK3oIJXV1ZPKHfnRKu9GGPiS0KCMHVidS3mheVb2HOkKHrvF7Wcm5g3N5RWLewzrHdbzurb3uOIjDHmeJMGZXJyV2c8enFZBY8v2lBLirqzAqYe8ovLmD5vLd/6/VyW7iyrOn7z2L4eRmWMMTUTEaZNqq7FPL98C1e/l8/p985l+ry15BeXBUkdHitg6ii/uIwLH1/KzMWbquYbA0gQ+N+3Vkf0h2SMMZE0tGcGqc2O/fg/WOBMznvh40sj9vllBUwdzVyykS0HCqq6JFeqUGeBn5lLott4ZowxdfXU+5uqHun7Ki6riOjnlxUwdfTS8i3HFS6VissqeGnF1hhHZIwxoXlp+RZKywOv4RjJzy8rYOootyD4XD65MRgla4wxdRGrzy8rYOooI+34NRWOPZ8co0iMMSY8sfr8sgKmjq4Y3pOUpMC3LyUpgSuG9YhxRMYYE5pYfX5ZAVNHU0b1oWe7tON+SClJCfRsl8aUUX1qSGmMMd6K1eeXFTB1lJ6SxOs3n82U0b1pm56MAG3Tk5kyujev33z2cSvGGWNMvIjV55d9CtZDekoS0yYOYNrEAXG7RKsxxgQSi88vq8EYY4yJCitgSVqh9QAACvtJREFUjDHGRIUVMMYYY6LCChhjjDFRIaqBpwtoikRkH7CljsnbA/sjGE6kWFzhsbjCY3GFpzHG1VNVA66uaAVMhIjISlUd4nUc/iyu8Fhc4bG4wtPU4rJHZMYYY6LCChhjjDFRYQVM5DzldQA1sLjCY3GFx+IKT5OKy9pgjDHGRIXVYIwxxkSFFTDGGGOiwgqYMIlIKxH5vYh8JCKHRWS3iLwuIv1DTD9IRBaISIGI7BSRe0UkMUKxTRaR2SKyS0RURK4OMd097vX+27e9jMtNe7aIfCgiRSKyWURui0RMPvnfICLr3fxXicj4ENJE7H7V9fdBRFqLyLMikuv+Hr4sIu3Cff9IxiUiWTXcl1kRjKuviMwUkS9EpFxEskNMF+37FXZc0b5fInKJiPxHRHaISJ77+/2jENKliMiDIrJXRPJF5G0RyapLDDabcvh6ADcATwP/A6QBvwI+FJFTVHVbTQlFJAOYD6wGzgf6AA/iFPS/jkBsFwNZwFvA9WGmPQz4f0B+E4GYoI5xiUhfYI6b7lfAGcB0ESlQ1b/XNyj3j+1J4B7gA+Aa4C0RGaqqX9WSvN73q56/D68C/XHuZwXwJ+ANYGQ4MUQhLoCfAUt99iM5sPBE4DxgBRB8WcZjRe1+1TMuiN79mgZsBqa6eZ4H/ENE2qvqI0HS/RXnb3YqsA/n72OeiJysqkVhRaCqtoWxAelAc79jbYE84He1pP0VkAu08jn2c6DA91g9YktwX1sAClwdYrp7gP1RvGd1jWsmsA5I8jn2OLANt4NKPeNaCzzjGyfwJfBSLO5XXX8fgOHufRzlc+wM99gED+PKcmP4brR/l9yv/wVkh5AmqverHnFF9X4B7QMc+wewOUiabkAZcKXPsa5ACXB9uDHYI7IwqWq+qhb6HTuIM8VMl1qSnwvMUdUjPsdmAc2B0RGIraK+eURDPeI6F5itqmU+x2bh/BGcVJ+YRKQ3zn+0r1Yec+N8zX3fWKjr78O5wB5VXVJ5QFU/wvlvNRKxR/X3tD7q+LsU7fsVl397qhqoJvQpwT+nJrmvs33y2YFTww/7XlkBEwEi0gHoi/PfdjADgTW+B1R1K85/hgOjE13I2ojIfhEpFZFPReQiL4MRkXSgO373i+rHUPW9X5XpA+Xf1v2ZBhOJ+1XX34fj0rm+qSVdtOOq9KzbDrFLRKaLSPMIxFQf0b5f9RXL+zWc4J9TA4Htqprnd7xO98raYCLjQZxHZM/Vcl0GcCjA8Vz3nFc24DwC+RRoCUwB/i0iP1DV2UFTRk8b99X/fuW6r/W9X5Xpg+W/r4a0kbpf/9/eucbaVVRx/PcXtFBTsIbIQxqKIkQ04YOPFMXQKqJERGshV201PrDRhBhjVBI1em3xgxA0xAeYijTGoMUHNkQv0mALbWg1QIJJY4kUFUtqDRZDUks1zfLDmls3557nPnvf29j/L5nss+fMzF5nnTln7T2zZlbd/tCv3stGuH7Tch0Cvg3cAzwDLAWuJedw3tWAXHVpW191mVV9FQeWdwMf6VOs0f8oGxjSwwQ4fVC5iJhxFyTpE8AqYEVE/ONokWsUIuKHHde9C3gA+BKVR+XZlmtUjlZ9HStExF7gmkrWFkn7gO9IuiAiHpkj0Y5KZlNfxQvsdmBjRKxvqt1B2MAkVwHrhiin55xIVwDfBK6NiDuHqP80cHKX/IX87855bLnGJSJC0s+Br0k6LiIOz4Fc03dRnfqavosaV1/T9U/muXds/drvyhD66sWo/aFar9sQ3qB6w1JXrm78lHTMeA0wVwambX01SeP6kvRiYIqcJ145oHiT373nYAAi4nsRoUGpWkfSG8mJz1si4oYhL7WLjnFMSYtIV+cZd9V15GqQKGnmG7MgV0QcIL3FOsd9e82djCrXro72qu3vj4hew2M9RaaHvvowUn/oV6/Qa65hVOrK1Y3oOM4FbeurSRrVl6T5pJv/C0hvtX8NqLILWFTmQKvU0pUNTA0kvQq4C7gbGGXh3xTwNkkLKnkTwEHgvuYkHA9JAlYAj4xwN94GU8DyjgV+E6ThGbROpS8R8Tg52XnVdJ6k55XzqVHaGkNfdfvDFHCapIsqMryWnE8YSfaG5erGleX4UANy1aVtfTVJY/qSdDzpFfkK4O0R8fchqt1Tjssr7ZxBrhcaXVej+jUf6wl4CfkH9wQ5Kbekks6vlDuLmf7kC4G9wCbgEmA16RxwXUOynU920FXkHdC3yvnFlTIXF7mqefeRhvLS0rF+RS5Gu2KO5Tqn6Od2YBk5sf4favjj95DrfcBhcvHgMtJJ4yDw6tnQ17D9gXQquLUj79fA48B7yInbR4GtDemlllzk+qAbi0yXAGuKPn/WhFzlGvNL37kS2A7srJzPnwt91ZWrbX2ROyRH6atLOtK8UuZe4N6Oet8lF2Z+gFxMvAP4I3DCyDI0peBjJZFGJXqkLZVyi+myqJD8s/1N6Uh7gbXAcQ3JNjmEXNPyL63k3Vp+fAeBA8BW4LIGdVZLrpJ/EfA74Fngz8AnG/4+P1Z++IeAh4G39Pi+W9HXMP2hfO71HXkvAm4j54+eIY3wjIV1Y+hlZLmA9wIPkrsc/LvodQ3lz6whuaZ/V93S4jnU18hyta2vcr1BMm2hY1EoMA/4OulFeYC8gTq7jgzert8YY0wreA7GGGNMK9jAGGOMaQUbGGOMMa1gA2OMMaYVbGCMMca0gg2MMcaYVrCBMaYGkpb2CHdbTR+aBTnWS3qw7esYUwdvdmlMPR4mY2t04xZyy/WtsyeOMUcfNjDG1CAy2uOOznxJq4ELgA9HxO5ZF8yYowgPkRnTEJLOA74BbIg+MTckTUr6W9lcs5r/jjK0dk45/6CkbZL2S3pa0uaySWM/GSYlzQiVW9q9piPvakk7JR2S9BdJnxvh4xozEBsYYxpA0vPJ/a2eAj4+oPgG4FRmxrefAB6KiMfK+WLgB+QOz+8nN1ndKmnsKIySPgvcDPwCuLy8XttphIwZBw+RGdMM15FDY0sjolvI2SNExB8k/Z40KJsBJM0jw+SurZRbM/26PO1sAl5P7kq9hppIOgn4Mrk78ldK9qYSO+SLkm6OuQ3TYP5P8BOMMWMiaRnwGeCrEbFtyGobgBUlZgfAZcAC4I5Ku6+UdGcJo3uYDFVwHnDumCJfCLwQ+Imk46cTuXvyqcCZY7ZvDGADY8xYSFpIDmP9ltGeKjYApwBvLucTwPaIeKK0u4AM/rQI+DQZ8Ol1ZBjdE8YU+5Ry3Ekarem0ueQvGrN9YwAPkRkzLuuAk4CVowwrRcTusn5lQtI24J3A5ytFLiSfJN4aEUdC1UrqFi+9yrNkeNwjFCNYZX85Xg7s69LGo4M/gTGDsYExpiaSPkqGSl4VEX+q0cSPgS+QQ1MnkuFtpzmxHA9VrvcGcuK/XzjdPcACSS+NiCdL3qUdZbaTgcTOiIhf1pDbmKGwgTGmBpJeDtxEroXZLWlJl2J7ImJPn2buAG4o6f6I2Ft5bwcZpnidpOvJp5lJ4MnORjq4mzQe35d0I3A2HV5tEfFPSZPATZLOAu4nh8vPBZZFxHKMaQDPwRhTjzeRE+VLyCeCbunqfg1ExF+BB4DTyaeZ6nv7SPfk04CNwKdIQ/EYfYiIp8inqjNJF+RVpItzZ7nrgdWkc8FG4EfASrz7gGkQh0w2xhjTCn6CMcYY0wo2MMYYY1rBBsYYY0wr2MAYY4xpBRsYY4wxrWADY4wxphVsYIwxxrSCDYwxxphW+C9ichEjGYkxTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot results for Z\n",
    "plt.plot(z_values, p_z, 'o-', linewidth=3, markersize=8)\n",
    "plt.grid()\n",
    "plt.xlabel('Z value', size=15)\n",
    "plt.ylabel('probability (%)', size=15)\n",
    "plt.title('Z Distribution', size=20)\n",
    "plt.xticks(size=15)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAElCAYAAADZb/T+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debgcVZ3/8feHJRL4QRIQA8gSFlmCcXQICBrgshuXQQNOhkUNipFxGFyiCAxCAEcBhwCKICgYYRgRJYhAkE1uEEFkFwxhk7AFQSAJhCwQ+P7+ONVJUenuW/emq2/u5fN6nn769qk6p071re5v1zl1TikiMDMza7WVersCZmbWPznAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgGmj5E0TFJImlxIn5ylD1uOsjuyMiZ2I8/ELE9HT7fbU5JmSprZ7u2WJWlVSSdIekTSoux9+mQbtx+SOtu1vSq14ziTNC7bxrhu5Fnmc9fKz2hPPpMrEgeYFsoOBA8sWkFlASlyj9clvSjpfkkXSfq0pAEt3OQE4DhgFvA/wAnAjBaW3209/aLOfWnmH4slPSfpakmjK6pyv9YoGPUXq/R2BaxljgZOBp5ZjjL+DGwDvNCSGq24zgTmkH5grQVsBXwKOBh4RNLBEfHnFmzn48A8YK+IeK0F5a0I5gJnZH+vBvwT8FHgo5K+EhE/6LWa9Y7ufO568hnt059JB5h+IiKeBZ5dzjLm08u/sNvkjIiYmU+QNAg4CfhP4DpJO0bE8r4XGwAv9qPgAjAnIibmEyQdAlwAfFfST7Pj6G2hO5+7nnxG+/pn0k1kFcufAmd/XyLpBUkLJd0p6eMN8q0paZKkp7N1Z0j6Og3+Z8X2XUk7Zq8vb1K3B7O+gbWz1w3beyVtJ+l3kl6R9LKkGyTt1NU+N1jeWWxKlDRA0uGSpkp6IqvXS9l2Km9+iYi5EXEEcCEwiPRLs1jv1SUdLeleSa9KmifpNkkHFNabnO3fpsAmuSalmbl1xkm6TNLfJC3I3tM/Sjq4Xv2a9TeVbfbK8h+fvbwp39zVLF8Jk4FXgTWAbbNt1Y7HzST9p6S/ZPvZmavPeyRdKOkZSa9JmpW9fk8X+/E5Sfdk5T0v6QJJ69VZbztJZ0q6LzuWFir1h50maUgX2/iYpFuz//NsSb+uV6/i566LMouf0YnA49niz+mtzY/jsnWafSbXlvS97HO8QNJcSTdK2rvOugMkHSHp7mx/5mfH1BWS9uyq7j3lM5j22YR0uvs34CJgbWAscIWkPSPiptqKkt4B3AhsD9wHXAwMBr4N7FpmYxHxJ0kPkZou1omIF/PLJe0AbA1cFhEvNStL0oeAG4ABwBTgUeD9QCfw+zL1KWFtUtPVrcD1wD+A9YFPAFMlfTEiftqibTVzIvBZ4OOS1oqIlwEkDSbt6weAu0m/2FcC9gH+T9K2EXFsVsZvgJnAV7PXtSalObntnAP8FbiZ9Kt2HVJT00WStoqIb1ewb2cAnyQdQz/P6tgqyp6LwepMYGfgamAq8AaApO1Jx9SawG+B6aTj8WBg3+wzcUed7XwN2Bv4JfA7YBRwCNAh6YMR8Y/cul8kNX1Oy7a1ErAd8HVgdLb+K3W2MQYYDVxOOsbfD+wH7CbpQxHxUJk3pIRO0uf6K6TP+W9yy+5tllHSJln+YcAfSO/FGqRm2d9J+lJE/CSXZTJwAPAA6UfUAtIZ9ijgI6T3p/Uiwo8WPUgfriikDaulA8cXlu2TpU8tpB+TpV8GrJRL3xR4KVs2uZBncpY+LJd2dJZ2eJ26/ihb9olcWkeWNjGXJtIpegD7Fsr4Sm7fOurs8+TidrPlnXXep3cAG9ZZdxDpQ/ESMLCwbCYwsxv/n5nF96jBek9l6+1W5/09srDuaqQP95vA+8vWD9i8TtoA0g+L14F3d6OsicX/Qe547Cyzbon3rvY/XaYOwOezZfNq/6Pc+/UMsGlhfQEPZssPKiwbm6XPKBz7tXq/BnygkOf0bNn5hfRNgJXr1PcL2frfKqSPyx3PH29wrN9Y4nNX9/jvzrrNPpO5z9CbwL8V0geTgtMCYGjuM/QmcGeD92Od7hwL3Xm4iax9ngC+k0+IiGuBJ4EdCuseQjogjoyIN3PrPw50pxP1oqycz+UTla6U+jfgeeCaLsr4EKkT/OaIuKKw7CzgsW7Up6GIWBQRT9dJn0s6WxhCOqNrh1on7LoAktYh/bK+MyJOLdRvIfAt0pfmgWU3EBHLvG+R+mp+RGpZ2KNHNa/e4KxJbqKkkyVNBc7Plh0TEQsK65+aHbd5HyKdrdwWERfnF0TEL4FbSMfcqDrbvygi7imkTSRdfHBgdvZfK+uJiHijThkXAC+TfuDV8/uIuKqQVjvWd8/OHnqNpH8inYVeFhGX5JdFxBxSM+hqpLMuSAFKwCLS9wGFPC8W01rFTWTtc2+Dg/0pYElfhqQ1gS2Ap+p9CZF+uRxfJ30ZEfG0pBuBvSQNj4jp2aJPkJqkTo+IxV0U88/Z87Q65b8h6RZg8zL16YqkbYFvAruQmsdWK6zy7lZsp0xVsudac8/2wMpAo/EIq2bP25TegLQxKTDtAWwMDCys0q597a5BLD3+3iCdWV4DnBURU+usX+9qvNox1ah59fek4PIBUhNiXr3jcK6ke0lfutuQNS9JWhX4EunH1PCs7vkf1Y3e466O9Q+QfjD2ltr3xaAGx+O62fM2ABHxsqQrSZ/7eyVdRmpWuz0qviDDAaZ95jRIX8xbD/pB2fNzDdb/eze3OxnYi3QW860srXZG8/MS+Vtdn7ok7Uj6YlmF1Ez0W9KvzDdJbeD7kprR2mGD7LnWnr9O9rw9zc+i/l+ZwiVtRvriHUL6oF9H+gX+BqnJ5HO0b1+764mIGNaN9esdH7VjqtEVVbX0wXWWdXUcDsql/ZLUB/M34IpsnUXZsq/S+D3uzjZ6Q+143Ct7NJI/HseSPv8HksZjASyU9GvgGxHRaJ+XiwPMimdu9jy0wfJlrpbpwuWkL+qDJR1DOjhHA/dFxH0V1ad2Gt7o+Kr3xXEs6Vf8bhHRmV8g6WhSgKmcpC2ADUmB/64sufYenB4RX2/BZr5O+j8cEhGTC9s/gEKTZuZNUh9NPfXezxVFvSvUau9no2N5/cJ6eV0dh3MBJI0kBZcbgNH5M3VJKwFHNqlzqW30otr2S487ypouJwITJW1EaiUYR2r6HUa6EKPl3Aezgol0VcujwLsl1Wt66uhmeQuAS0m/yvck/YJZhXJnL5CumII6V69JWpn67eSzs+eN6uRZC9iyTp4tgJeKwaXRtit0XPZ8ZSy9wujPpC/4Vn0It8ieL6uzrNG+zgaGZs0+RSO7se1aM+3K3cjTarU+lI4Gy3fLnu+us6zecTiIdJa7kHTxACx9j39bpxl4B5ZtkuxqG/ljvdgHtDx68v/4U/bco+MxIp7K+r72IX3XjMr6GVvOAWbF9DPS/+aU7NcWAJI2BY7oQXmTs+fPZo/FpEufy7gVeAjYRVLxLOJw6vS/ZF/MM4APSxpeS88+pJOo/+GeCawt6X35RElfoHFnbMtIWkvSD4DPkJozj6oti4jnSe/XSEnfzvajmH/z7P9TxszsuaNQxj7AoQ3y/Jn0w+CQQp5xwIdLbheg1qG7cTfytNofScfUKEn75xdkr3cGHiZ19hd9RtIHCmkTSc1Wv4iIWhPYzOy5o1D+u0gXUjSzu5Ydn1Y71m+KiFb2v8wmneWV/n9ExJ2kptUxkj5fbx1JI7J9RdK6kkbUWW0NUjPaYtLVeS3nJrIV02mk8Qr7AXdLupbUDPKvpE7Pf+lOYRHxR0mPAp8mdUhfmX1plskb2Zf89cBlkvLjYPYgXaL7kTpZv0+6uuiPkn5F+nW5W7b9+0hTjOSdQQokt0i6lNQMMJL0q/HXwP60zlclzSF15temitmF9IF7GDg4Ih4u5DkceA9pnMxnsg7f50hnhtuQ+mYOYOnAuWbOJgWKX2Vt4LOA95Lex0tJ7eVFP8zynCNpD9LFIe8ndfheRRr/UMZNpLOx70l6L9nZZkR8p2muFsqOqc+RjqlfSrqC9INkK9Jx/wrw2fwVlDnXkI6pS0l9NaOyx0xyPwqAO0iBbIykW0nBaiipefgh0nveyJXA5UqDlGvH+mjSBQ1f7sk+NxIR8yTdDuws6WLS8fcG6czrL02yHkjqszxf0hHA7aQfRhsC7yMdTzuRrhR9N3CPpPuBv5COnbVIx8x6wA+i/nig5VfV9c9vxwfNx8FMbpCns5gnS1+L9Gv/GdKX8wzS5Imb1SuPOtfYF5YfW6sfsF+DdTqoc819tmw7UjB5JXvcQDqAJ9JgXAVpvMFfSR2rfwfOJfU9NNrnj5NO/18hfViuY2lbcQDjCuvPpGfjYGqP10lfGveTLuneHxjQJP8AUqC5lRQAF5EuM7+R1Gm8Ttn6kS7V/T3pC/4V0hfgJ7v4H4wi/cCYT+pXu5r0ZVL3f0CdcTBZ+sEsHSuxzDHboL6147jU+93V8Zits1X2vj+b/S+eBf4X2KrOukv2MTseavX/B+mMf/06edYmBfOZpM/QY8B3gdXr/W/yx1l2LN5GmqFgDqk5c8sy+0k3xsFk6VuQgtqLpOC/5Fjv4nhYkzRm7i7SGKQFpB84VwPjgTWy9QaTmn5/T/o+WZS9152kH0Uq+xnq7kNZBczMzFrKfTBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhy5Rz3vnOd8awYcN6uxp93quvvsoaa6zR29Uwa8rHaWvcddddL0TEuvWWOcDkDBs2jDvvvLO3q9HndXZ20tHR0dvVMGvKx2lrSGo48NRNZGZmVgkHGDMzq0TbA4yk4dl9o+cr3YP7xHpzOxXybC/pZ5IezfI9JOl4SasV1qvdm7z4qDeViZmZVaitfTCShpCmGJlOmn59c9K8WyuRpjJpZGy27inAI6TpMU7KnvcrrDuXZefGehAzM2urdnfyH0aaSXdMRLwMXJ9N3z5R0qlZWj0nR8QLudedkhYC50raJN46u+niiPgTZmbWq9rdRDYauLYQSC4hBZ2G9/woBJea2j0ZNqizzMzMelm7A8zWpFmBl4iIJ0mzw27dzbJ2Is08Wrxv/WBJL0h6XdI9ksb0uLZmZtZj7W4iG0L9e9PPzpaVImk9Up/NRfHW+5o8SroV6j2kqay/RLqHyX4RMaVBWeNJU1szdOhQOjs7y1bDGpg3b57fR1vh+TitXlun65f0OvDNiDijkP40cGFEHFOijAGkCwU2BLaLiNlN1hXp3h0DI+L9XZU9cuTI8EDL5ecBbNYX+DhtDUl3RUTd23a3+wxmNunWpkVDWHof94aygHEhsC3w4WbBBZbcOW8K6dbDK0fEG83WN+vPhh11dW9XYYUyYcRixvk9AWDmyR+rpNx2B5gZFPpaJG1EusPcjLo53uoM0uXNe0VEmfVh6d0LzcysjdrdyX8NsI+kNXNpY0m3+pzWLKOko0m3qz04Im4ps7HsjGc/4D6fvZiZtVe7z2B+DBwBTJF0Cun+8hOBSflLlyU9CkyLiC9krw8k3Ut7MvCMpB1zZT4WEf/I1ptGunf2DGAN4IvAB0n3OjczszZqa4CJiNmS9gDOAq4kXVF2OinIFOuVnz5m7+x5XPbIO4QUeCBdRfZVYH3SJcx3Ax+LiGtaUX8zMyuv7dP1R8R0YPcu1hlWeD2OZQNLvXxfWI6qmZlZC3k2ZTMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSrQ9wEgaLulGSfMlzZJ0oqSVu8izvaSfSXo0y/eQpOMlrVZn3Q9Lul3SQkmPSzqiur0xM7NGVmnnxiQNAW4ApgP7ApsDp5EC3bFNso7N1j0FeAR4H3BS9rxfrvwtgGuBq4CjgR2ASZLmR8RPW70/ZmbWWFsDDHAYMBAYExEvA9dLWguYKOnULK2ekyPihdzrTkkLgXMlbRIRT2Tp3wRmAQdHxGLg95I2Bo6XdH5ERDW7ZWZmRe1uIhsNXFsIJJeQgs6ujTIVgkvNPdnzBoXyp2TBJV/+hsB7e1RjMzPrkXYHmK2BGfmEiHgSmJ8t646dgDeBxwAkrQFsVCwfeDC3bTMza5N2N5ENAebUSZ+dLStF0nqkPpuLIuL5LHlw9lwsf3Zu2/XKGg+MBxg6dCidnZ1lq2ENzJs3z+/jCmjCiMVdr/Q2MnSg35Oaqj6v7Q4wy03SAOBSYB7wteUtLyLOA84DGDlyZHR0dCxvkW97nZ2d+H1c8Yw76urersIKZcKIxZx2f5/7CqzEzIM6Kim33e/ubGBQnfQhLD3TaEiSgAuBbYEPR0Q+T+3MpVh+7cyly/LNzKx12h1gZlDoC5G0EbA6y/ad1HMG6fLmvSKi2JfzqqSniuXnXpcp38zMWqTdnfzXAPtIWjOXNhZYAExrllHS0cDhpEuQb2lS/qcKAzfHAk8BD/S41mZm1m3dCjCSBkh6t6TNs0GT3fVjYBEwRdKeWQf7RGBS/tLlbMT++bnXBwLfJTWPPSNpx9xj3Vz53yddknyRpN0kHQl8CTjRY2DMzNqryyYySdsCnwX2BEYAK+eWvQjcCvwauCwiFjQrKyJmS9oDOAu4ktRvcjopyBTrlT8L2Tt7Hpc98g4BJmflPyrpI8Ak0tnM34EJHsVvZtZ+DQOMpA8D3wF2Ae4gNWH9AHiBdBYyGBgGjCQFiR9KmgScHhHzGpUbEdOB3ZtVKiKGFV6PY9nA0ijvLaQpYszMrBc1O4OZQgoon4mIp5sVkvV57Al8NUs6qTXVMzOzvqpZgNkkIhaWKSQi3iBNMnltvRmOzczs7adhJ3/Z4NKqfGZm1r90exyMpFVIU6vsBgi4CTi3MMGkmZm9zfVkoOWZwHbARcAawJGk+7J8qYX1MjOzPq7ZVWT/HBF311m0P/Ce2rgVSX8lBRsHGDMzW6LZQMtrJP2kMJAR0tiSvWDJ3GC7Ac9WVD8zM+ujmgWYrUj3aXlQ0jclrZqlH04a8/I8aQLJA4F/r7aaZmbW1zS7imxORHwF2Jl0ljJd0r9ExB9IAyz3JA3C3DQibm5HZc3MrO/oci6yiHgwIj5KGkR5iqTrgS0i4i/ZY1HltTQzsz6nywAjaaCkQRFxNWkusmuAaZLOkrR25TU0M7M+qWGAkbSFpD8CrwIvSXoE2DUiJgHDSVegzZD0n4Xp8c3MzJqewVwIzATWI01s+XPSNPurRcQ/IuIw0tVknwLur7qiZmbWtzQLMMOByRHxfES8ApwNrEm63woAEXFfROwO/Fe11TQzs76m2Uj+60id+muR7jg5Hngse7xFRFxeTfXMzKyvahZgPg8cCxwFDADuAvb0nSHNzKyMhgEmu2nYUW2si5mZ9SPNriJTTwrsaT4zM+tfmnXyPyzpUElrlClI0naSLsRnPWZmRvM+mG8BJwBnSroOuBV4AHgBWES6dHlT0tT9HwE2An4KXFBlhc3MrG9o1gczRdLlpDnHPkuaKmZ9oNbJL+A1Uuf/ucBFEfF8tdU1M7O+oukNx7Irxq7PHkhanzTwcjXgJWCm5yIzM7N6unVHy4h4Ft/7xczMSuhyskszM7OecIAxM7NKOMCYmVklHGDMzKwSpQKMpBFVV8TMzPqXsmcw90m6Q9K/Sxq8PBuUNFzSjZLmS5ol6cSublgmaYCk70v6g6QFkupOuClpsqSo89h6eepsZmbdVzbA7A5MB04FZkn6haS9ujvvmKQhwA2kwZr7AicCE0gzBjSzOnAoMJ80o0AzM4CdCo+Z3amnmZktv1LjYCKiE+iU9GVgLDAOuBZ4WtLPSTcmW+Y+MXUcBgwExkTEy8D12f1mJko6NUurt/05ktaOiJB0OCngNfJqRPypzH6ZmVl1utXJHxGvRsQFEbELsBXpzOAY0sSY0yR9qosiRgPXFgLJJaSgs2sX2/Z9aMzM+pBuX0UmaZikiaQzmJ2AqaS7XT4H/FLS6U2yb01qwloiIp4kNX21qp9kuKSXJS2SdIukpoHLzMyqUaqJTNLqwP7AIcDOwOPAT0hNY7WpY86XdAhwJvC1BkUNAebUSZ+dLVte9wC3k/qL1iX171wvaVRE/LleBknjSQGSoUOH0tnZ2YJqvL3NmzfP7+MKaMKIxb1dhRXK0IF+T2qq+ryWnYvsOdLZzhTSbZMb1eYO4MUW1KtHIuLM/GtJU4G/kprxPtkgz3nAeQAjR46Mjo6OimvZ/3V2duL3ccUz7qire7sKK5QJIxZz2v3dmo6x35p5UEcl5ZZ9d48E/i8i5jZbKSIeIN0jppHZwKA66UOyZS0VEfOzIPOJVpdtZmbNle2DWReoe2dLSetLOq5kOTMo9LVI2oh0GfKMujmWX7D0HjZmZtYmZQPM8cCGDZZtkC0v4xpgH0lr5tLGAguAaSXLKE3SQOBjpJuimZlZG5VtIhONzwI2pHzz1o+BI4Apkk4BNgMmApPyly5LehSYFhFfyKWNJp1FvT97vX+26I6IeELSIOAq4H+BR4F3ki422AD4dMn69dgwt28vMWHEYrf358w8+WO9XQWzXtEwwEj6HPC57GUA50gqDoRcDRgBXFdmYxExW9IewFnAlaQryk4nBZlivYrTx5wDbJJ7/avs+RBgMrAI+AdwLPAuYCFwG7BrRNxZpn5mZtY6zc5g5rP0ijABc0m3Sc57jdTsdXbZDUbEdJqPxCcihpVJKyxfCIwpWw8zM6tWwwATEb8iO0uQ9DPgpIj4W7sqZmZmfVvZucgOqboiZmbWvzTrgzkV+EFEPJ393UxExLdaWzUzM+vLmp3BfBq4GHiarq/CCsABxszMlmjWB7Npvb/NzMzK6PZsymZmZmU064P5aHcKioipy18dMzPrL5r1wVxF6lspc1vkYNmBkWZm9jbWLMC438XMzHqsWSf/E+2siJmZ9S/N+mBWj4j5tb+7Kqi2rpmZGTRvIntF0k7ZrYbn0fU9VdwHY2ZmSzQLMJ8HHsv97Zt2mZlZac36YH6e+3tyW2pjZmb9RtkbjgEgaTDwXmB94FnggYiYU0XFzMysbysVYCStAvw38B9AvsN/vqSzgf+KiNcrqJ+ZmfVRZc9gJgHjgROBKcDzpLtG7ke6g+RqpFshm5mZAeUDzGeAYyJiUi7tJeC/JS0kBRkHGDMzW6LsZJdvAn9tsOwBfIWZmZkVlA0wFwGHNlj2ReB/W1MdMzPrL5qN5P9y7uVMYH9JfwV+y9I+mH2BNYH/qbCOZmbWBzXrgzmrTtoGwDZ10icBZ7akRmZm1i80G2jpm5GZmVmPOYiYmVklujuSf0NgS9K4l7fwHS3NzCyv7Ej+NYFLgb1rSdlz/vJkz6ZsZmZLlG0i+x6wMbAzKbh8CugAzgceB3asonJmZtZ3lQ0wHyXNRXZ79npWRNwcEeOBK4BvVlE5MzPru8oGmKHAUxHxBvAqsHZu2VSWNp2ZmZkB5QPMU8A7s78fAT6eW/ZBYGHZDUoaLulGSfMlzZJ0oqSm/TeSBkj6vqQ/SFogqeHUNJL2lXS/pIWSpksaW7ZuZmbWOmUDzPXAntnfpwP/IelWSTcBJwEXlilE0hDgBtLFAfuSZmeeAJzQRdbVSVPVzAdubVL+KOAy4CZgNHA18AtJPsMyM2uzspcpf4vsPjARcZGkecD+wEDgcODckuUcluUZExEvA9dLWguYKOnULG0ZETFH0toREZIOB3ZvUP63gZsjojaz802StgWOA64rWUczM2uBUmcwETE/Il7Ivb48Ig6KiDERcU5EvFlye6OBawuB5BJS0Nm1izo0nbFZ0juA3UiXU+ddAuwkaVDJOpqZWQt0ayS/pK0kHSzpm5IOkrRVN7e3NTAjnxART5KavrbuZllFmwOrFssHHiTt55bLWb6ZmXVD2YGWawE/Id3BciVgHvD/gDclTQEObdS8VTAEmFMnfXa2bHnU8hfLn11Y/haSxpPu1snQoUPp7Ozs0cYnjFjco3z90dCBfj/yenpMtZr/J2/l43Spqo7Rsn0wZ5MuRf4scHlELJA0EBhDmnX5bODgSmpYsYg4DzgPYOTIkdHR0dGjcsYddXULa9W3TRixmNPu79YsRP3azIM6ersKgI/RIh+nS1V1jJZ9d/cFvhYR/1dLiIgFwMWSVidN11/GbKBeX8gQlp5p9FQtf7H8IYXlZmbWBmX7YOYBzzZYNos0+LKMGRT6WiRtRLpCrdh30l2PAa8Xy89evwk8vJzlm5lZN5QNMD8CvpE1iy2Rnb18g9REVsY1wD7Z5Jk1Y4EFwLSSZdQVEYtI418+XVg0FrgtIuYuT/lmZtY9zW6ZfGoh6T3AU5KuZ+ktk/ciBYc7S27vx8ARwBRJpwCbAROBSfmLBCQ9CkyLiC/k0kYDawDvz17vny26IyKeyP4+CeiUdAbwG9Icah8FPlKyfmZm1iLN+mCKZwKvZ4/8zMmvZM/7UWLCy4iYLWkP0oUBV5Ku+DqdFGSK9SpOH3MOsEnu9a+y50OAyVn5t2SB5zvAv5Nmej4wIjzI0syszZrdMnnTKjYYEdNpPBK/ts6wMmkN8v6GdPZiZma9yLdMNjOzSpQOMJI2k3RONlPxM9nz2ZI2q7KCZmbWN5Udyb8d6QqthcBVwHOke8TsBxwkabeIuLuyWpqZWZ9TdqDl/wD3AKMjYn4tMbtMeWq2vGm/ipmZvb2UbSLbATg1H1wgzbJMCi4fbHXFzMysbysbYBYA6zRYtjbduKOlmZm9PZQNMFcDJ2d3jFwie/090pgWMzOzJcr2wXwduAKYJul5lo7kfxdwG+m2x2ZmZkuUCjAR8SIwStJHgO2B9UmTX97uUfJmZlZPlwFG0mrAX4AjIuJ3wO8qr5WZmfV5XfbBRMRCYDBpynszM7NSynbyX0yaVNLMzKyUsp38TwL/KukO0j1dngMitzwi4pxWV87MzPqusgHmtOx5fWC7OrNpAzoAAAqnSURBVMuDNJ2+mZkZUP4qMs+6bGZm3eLAYWZmlSjbRIakAcA40rxkS8bBAD+PiNcqqZ2ZmfVZpc5gJG0DPAL8CHgv8Eb2/CPgUUnDK6uhmZn1SWXPYM4D5gI7R8STtURJG5PuD/NjYJfWV8/MzPqqsn0wI4Hj8sEFIHt9PGn6GDMzsyXKBpiZwGoNlq1GGidjZma2RNkAcxTwHUlvubGYpB2Bk4BvtbpiZmbWt5XtgzkWWAu4tc50/S8Cx0g6prZyROzQ6oqamVnfUjbAPJA9zMzMSik7kt8TXZqZWbd4JL+ZmVXCAcbMzCrhAGNmZpVwgDEzs0q0PcBIGi7pRknzJc2SdKKklUvkGyTpZ5JmS5or6WJJ6xTWmSwp6jy2rm6PzMysntKzKbeCpCHADcB0YF9gc9LNzFYijbVp5lJgS+BQ4E3gFOA3wM6F9Waw7O2dZy5Pvc3MrPvaGmCAw4CBwJiIeBm4XtJawERJp2Zpy5C0E7A3sGtE3JylPQPcLmnPiLght/qrEfGnanfDzMy60u4mstHAtYVAcgkp6OzaRb7nasEFICL+DDyeLTMzsxVMuwPM1qQmrCWyGZnnZ8tK58s8WCffcEkvS1ok6RZJzQKXmZlVpN1NZEOAOXXSZ2fLepJvs9zre0h32ZwOrAtMIDXDjcrOeJYhaTwwHmDo0KF0dnZ2sQv1TRixuEf5+qOhA/1+5PX0mGo1/0/eysfpUlUdo+0OMJWKiDPzryVNBf4KHAN8skGe80g3VGPkyJHR0dHRo22PO+rqHuXrjyaMWMxp9/erQ2u5zDyoo7erAPgYLfJxulRVx2i7m8hmA4PqpA/JlrU0X0TMB6YC/9yNOpqZWQu0O8DMoNBnImkjYHXq97E0zJdp1DeTF9nDzMzaqN0B5hpgH0lr5tLGAguAaV3kW0/SqFqCpJGk/pdrGmWSNBD4GHDX8lTazMy6r90B5sfAImCKpD2zDvaJwKT8pcuSHpV0fu11RNwGXAdcKGmMpE8CFwO31MbAZCP9/yDpS5L2kDQWuAnYAPhuu3bQzMyStvZwRcRsSXsAZwFXkq4MO50UZIr1Kk4fMzZb9wJSYLwKOCK3fBHwD9KMAO8CFgK3kQZn3tnSHTEzsy61/RKKiJgO7N7FOsPqpM0hTQFT9+ZnEbEQGNOCKpqZWQt4NmUzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0o4wJiZWSUcYMzMrBIOMGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYMzOrhAOMmZlVwgHGzMwq4QBjZmaVcIAxM7NKOMCYmVklHGDMzKwSDjBmZlYJBxgzM6uEA4yZmVXCAcbMzCrhAGNmZpVwgDEzs0q0PcBIGi7pRknzJc2SdKKklUvkGyTpZ5JmS5or6WJJ69RZb19J90taKGm6pLHV7ImZmTXT1gAjaQhwAxDAvsCJwATghBLZLwU6gEOBccD2wG8K5Y8CLgNuAkYDVwO/kLR3S3bAzMxKW6XN2zsMGAiMiYiXgeslrQVMlHRqlrYMSTsBewO7RsTNWdozwO2S9oyIG7JVvw3cHBFHZK9vkrQtcBxwXXW7ZWZmRe1uIhsNXFsIJJeQgs6uXeR7rhZcACLiz8Dj2TIkvQPYjXSmk3cJsJOkQctffTMzK6vdAWZrYEY+ISKeBOZny0rnyzyYy7c5sGqd9R4k7eeWPaivmZn1ULubyIYAc+qkz86W9STfZrl1qLPe7MLyt5A0HhifvZwn6aEm9bASjoB3Ai/0dj1WFDqlt2tg9fg4XWo5j9FNGi1od4BZ4UTEecB5vV2P/kTSnRExsrfrYdaMj9PqtbuJbDZQry9kCEvPNHqar/ZcXG9IYbmZmbVBuwPMDAp9LZI2Alanfh9Lw3yZfN/MY8DrddbbGngTeLgH9TUzsx5qd4C5BthH0pq5tLHAAmBaF/nWy8a5ACBpJKn/5RqAiFhEGv/y6ULescBtETF3+atvJbnJ0foCH6cVU0S0b2NpoOV04AHgFFKAmAScERHH5tZ7FJgWEV/IpV0LvAf4BumM5BTg+YjYObfOKKATOIs0CPOj2fofiQiPgzEza6O2nsFExGxgD2Bl4ErSCP7TgeMLq66SrZM3lnSWcwFwIXAX8KlC+bcA+wN7AtcC/wIc6OBiZtZ+bT2DMTOztw/Ppmwt09OJTM3aQdIWks6V9BdJb0jq7O069Xdv+3Ew1hq5iUynkyYy3Rw4jfQj5tgmWc3aZVtSv+yfSLN+WMXcRGYtIelo4Ehgk9pcc5KOBCYC6zWayNSsXSStFBFvZn//GnhnRHT0bq36NzeRWav0dCJTs7aoBRdrHwcYa5WeTmRqZv2UA4y1Sk8nMjWzfsoBxszMKuEAY63S04lMzayfcoCxVunpRKZm1k85wFir9HQiUzPrpzzQ0lrlx8ARwBRJtYlMJwKTPAbGVgSSVicNtAR4N7CWpP2z11MjYn7v1Kz/8kBLaxlJw0kzWe9EuqLsp8DEiHijVytmBkgaBjzeYPGmETGzbZV5m3CAMTOzSrgPxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEg4wZmZWCQcYszaQ9LikkLRFL2x7B0kT271dMwcYs4pJ2gkYlr08oBeqsANwfC9s197mHGDMqncA8CpwO70TYMx6hQOMWYUkrQz8K/Bb4AJgG0n/lFs+WNJPJc2StFDSk5J+klu+oaRLJT0vaYGkxySdVNjGzpKmSZov6UVJP6lNOippHPDD7O/IHp2V77gZnuzSrGq7AUOBS4BbSHO1HQDcly2fBHwI+Brwd2AjYJdc/guBgcB40vxum5G7LYKkDwM3AL8B9gfWAU4m3Ydnf+Bq4DRgAmmOOABPPmpt4bnIzCok6XxgDDA0Il6TdBXwXtLkiiHpAeDciPhhg/zzgAMi4soGy/8ALI6I3XJpuwM3AiMi4gFJhwM/jAi1du/MmnMTmVlFJA0gBZfLI+K1LPkSYBOWnk3cC3xT0pclbVmnmHuB70kaJ2njQvmrZ+VcKmmV2oN0pvQ6sF3r98qsPAcYs+qMBgYDU7O+lsFAJ7CIpZ39h5Oat44DHpL0iKR/y5UxFrgTOB14QtK9kvbIlg0BVgbOJgWU2mMRsCqpuc2s17iJzKwiki4hBYh6ngPenb9XjqT3AUeSgs+IiJieW7YS6XLjicDOwMbAQuCVLG1qnW3MiohZbiKz3uIAY1YBSWsAzwNXAOcVFn+A1Lm/d0RcX8i3PjAL2C8iptQpdyfgVmC7iLhb0q3A4xFxUJO6jAfOBQZGxMLl2C2zbnGAMauApAOBi4EdI+L2wrJVgWdJly5vCVwOPAAE8EVS09rWpLOTa0lXkj0MvIN0Ndg2wGYRsUDSKFKH/qXAr7M8GwMfA/4rIh6WtAswDTgK+D3wckQ8VN3emyUOMGYVkHQlsFVE1Ou4R9LZwIHAz4A9SSP93wDuAY6LiD9IegfpsuZdSP0p84E/AUdHxP25sj4InEC63Hll4Angd8AJETFXkoBTgIOB9YCbI6Kj1ftsVuQAY2ZmlfBVZGZmVgkHGDMzq4QDjJmZVcIBxszMKuEAY2ZmlXCAMTOzSjjAmJlZJRxgzMysEv8fNeptZaNQCFAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot results for default probabilities\n",
    "plt.bar(range(K), p_default)\n",
    "plt.xlabel('Asset', size=15)\n",
    "plt.ylabel('probability (%)', size=15)\n",
    "plt.title('Individual Default Probabilities', size=20)\n",
    "plt.xticks(range(K), size=15)\n",
    "plt.yticks(size=15)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expected Loss\n",
    "\n",
    "To estimate the expected loss, we first apply a weighted sum operator to sum up individual losses to total loss:\n",
    "\n",
    "$$ \\mathcal{S}: |x_1, ..., x_K \\rangle_K |0\\rangle_{n_S} \\mapsto |x_1, ..., x_K \\rangle_K |\\lambda_1x_1 + ... + \\lambda_K x_K\\rangle_{n_S}. $$\n",
    "\n",
    "The required number of qubits to represent the result is given by\n",
    "\n",
    "$$ n_s = \\lfloor \\log_2( \\lambda_1 + ... + \\lambda_K ) \\rfloor + 1. $$\n",
    "\n",
    "Once we have the total loss distribution in a quantum register, we can use the techniques described in [Woerner2019] to map a total loss $L \\in \\{0, ..., 2^{n_s}-1\\}$ to the amplitude of an objective qubit by an operator\n",
    "\n",
    "$$ | L \\rangle_{n_s}|0\\rangle \\mapsto \n",
    "| L \\rangle_{n_s} \\left( \\sqrt{1 - L/(2^{n_s}-1)}|0\\rangle + \\sqrt{L/(2^{n_s}-1)}|1\\rangle \\right), $$\n",
    "\n",
    "which allows to run amplitude estimation to evaluate the expected loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add Z qubits with weight/loss 0\n",
    "from qiskit.circuit.library import WeightedAdder\n",
    "agg = WeightedAdder(n_z + K, [0]*n_z + lgd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.circuit.library import LinearAmplitudeFunction\n",
    "\n",
    "# define linear objective function\n",
    "breakpoints = [0]\n",
    "slopes = [1]\n",
    "offsets = [0]\n",
    "f_min = 0\n",
    "f_max = sum(lgd)\n",
    "c_approx = 0.25\n",
    "\n",
    "objective = LinearAmplitudeFunction(\n",
    "    agg.num_sum_qubits,\n",
    "    slope=slopes, \n",
    "    offset=offsets, \n",
    "    # max value that can be reached by the qubit register (will not always be reached)\n",
    "    domain=(0, 2**agg.num_sum_qubits-1),  \n",
    "    image=(f_min, f_max),\n",
    "    rescaling_factor=c_approx,\n",
    "    breakpoints=breakpoints\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the state preparation circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">             ┌───────┐┌────────┐      ┌───────────┐\n",
       "    state_0: ┤0      ├┤0       ├──────┤0          ├\n",
       "             │       ││        │      │           │\n",
       "    state_1: ┤1      ├┤1       ├──────┤1          ├\n",
       "             │  P(X) ││        │      │           │\n",
       "    state_2: ┤2      ├┤2       ├──────┤2          ├\n",
       "             │       ││        │      │           │\n",
       "    state_3: ┤3      ├┤3       ├──────┤3          ├\n",
       "             └───────┘│  adder │┌────┐│  adder_dg │\n",
       "objective_0: ─────────┤        ├┤2   ├┤           ├\n",
       "                      │        ││    ││           │\n",
       "      sum_0: ─────────┤4       ├┤0 F ├┤4          ├\n",
       "                      │        ││    ││           │\n",
       "      sum_1: ─────────┤5       ├┤1   ├┤5          ├\n",
       "                      │        │└────┘│           │\n",
       "    carry_0: ─────────┤6       ├──────┤6          ├\n",
       "                      └────────┘      └───────────┘</pre>"
      ],
      "text/plain": [
       "             ┌───────┐┌────────┐      ┌───────────┐\n",
       "    state_0: ┤0      ├┤0       ├──────┤0          ├\n",
       "             │       ││        │      │           │\n",
       "    state_1: ┤1      ├┤1       ├──────┤1          ├\n",
       "             │  P(X) ││        │      │           │\n",
       "    state_2: ┤2      ├┤2       ├──────┤2          ├\n",
       "             │       ││        │      │           │\n",
       "    state_3: ┤3      ├┤3       ├──────┤3          ├\n",
       "             └───────┘│  adder │┌────┐│  adder_dg │\n",
       "objective_0: ─────────┤        ├┤2   ├┤           ├\n",
       "                      │        ││    ││           │\n",
       "      sum_0: ─────────┤4       ├┤0 F ├┤4          ├\n",
       "                      │        ││    ││           │\n",
       "      sum_1: ─────────┤5       ├┤1   ├┤5          ├\n",
       "                      │        │└────┘│           │\n",
       "    carry_0: ─────────┤6       ├──────┤6          ├\n",
       "                      └────────┘      └───────────┘"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define the registers for convenience and readability\n",
    "qr_state = QuantumRegister(u.num_qubits, 'state')\n",
    "qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')\n",
    "qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')\n",
    "qr_obj = QuantumRegister(1, 'objective')\n",
    "\n",
    "# define the circuit\n",
    "state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')\n",
    "\n",
    "# load the random variable\n",
    "state_preparation.append(u.to_gate(), qr_state)\n",
    "\n",
    "# aggregate\n",
    "state_preparation.append(agg.to_gate(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
    "\n",
    "# linear objective function\n",
    "state_preparation.append(objective.to_gate(), qr_sum[:] + qr_obj[:])\n",
    "\n",
    "# uncompute aggregation\n",
    "state_preparation.append(agg.to_gate().inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
    "\n",
    "# draw the circuit\n",
    "state_preparation.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we use QAE to estimate the expected loss, we validate the quantum circuit representing the objective function by just simulating it directly and analyzing the probability of the objective qubit being in the $|1\\rangle$ state, i.e., the value QAE will eventually approximate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact Expected Loss:   0.6409\n",
      "Exact Operator Value:  0.3906\n",
      "Mapped Operator value: 0.6640\n"
     ]
    }
   ],
   "source": [
    "# evaluate resulting statevector\n",
    "value = 0\n",
    "for i, a in enumerate(job.result().get_statevector()):\n",
    "    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]\n",
    "    am = np.round(np.real(a), decimals=4)\n",
    "    if np.abs(am) > 1e-6 and b[0] == '1':\n",
    "        value += am**2\n",
    "\n",
    "print('Exact Expected Loss:   %.4f' % expected_loss) \n",
    "print('Exact Operator Value:  %.4f' % value)\n",
    "print('Mapped Operator value: %.4f' % objective.post_processing(value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we run QAE to estimate the expected loss with a quadratic speed-up over classical Monte Carlo simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:    \t0.6409\n",
      "Estimated value:\t0.6913\n",
      "Confidence interval: \t[0.6224, 0.7601]\n"
     ]
    }
   ],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae = IterativeAmplitudeEstimation(state_preparation=state_preparation,\n",
    "                                  epsilon=epsilon, alpha=alpha,\n",
    "                                  objective_qubits=[len(qr_state)],\n",
    "                                  post_processing=objective.post_processing)\n",
    "result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)\n",
    "\n",
    "# print results\n",
    "conf_int = np.array(result['confidence_interval'])\n",
    "print('Exact value:    \\t%.4f' % expected_loss)\n",
    "print('Estimated value:\\t%.4f' % result['estimation'])\n",
    "print('Confidence interval: \\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cumulative Distribution Function\n",
    "\n",
    "Instead of the expected loss (which could also be estimated efficiently using classical techniques) we now estimate the cumulative distribution function (CDF) of the loss.\n",
    "Classically, this either involves evaluating all the possible combinations of defaulting assets, or many classical samples in a Monte Carlo simulation. Algorithms based on QAE have the potential to significantly speed up this analysis in the future.\n",
    "\n",
    "To estimate the CDF, i.e., the probability $\\mathbb{P}[L \\leq x]$, we again apply $\\mathcal{S}$ to compute the total loss, and then apply a comparator that for a given value $x$ acts as\n",
    "\n",
    "$$ \\mathcal{C}: |L\\rangle_n|0> \\mapsto \n",
    "\\begin{cases} \n",
    "|L\\rangle_n|1> & \\text{if}\\quad L \\leq x \\\\\n",
    "|L\\rangle_n|0> & \\text{if}\\quad L > x.\n",
    "\\end{cases} $$\n",
    "\n",
    "The resulting quantum state can be written as\n",
    "\n",
    "$$ \\sum_{L = 0}^{x} \\sqrt{p_{L}}|L\\rangle_{n_s}|1\\rangle + \n",
    "\\sum_{L = x+1}^{2^{n_s}-1} \\sqrt{p_{L}}|L\\rangle_{n_s}|1\\rangle, $$\n",
    "\n",
    "where we directly assume the summed up loss values and corresponding probabilities instead of presenting the details of the uncertainty model.\n",
    "\n",
    "The CDF($x$) equals the probability of measuring $|1\\rangle$ in the objective qubit and QAE can be directly used to estimate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">                               \n",
       "  state_0: ──■──────────────■──\n",
       "             │              │  \n",
       "  state_1: ──┼────■─────────┼──\n",
       "             │  ┌─┴─┐┌───┐  │  \n",
       "compare_0: ──┼──┤ X ├┤ X ├──┼──\n",
       "           ┌─┴─┐└─┬─┘└───┘┌─┴─┐\n",
       "     a0_0: ┤ X ├──■───────┤ X ├\n",
       "           └───┘          └───┘</pre>"
      ],
      "text/plain": [
       "                               \n",
       "  state_0: ──■──────────────■──\n",
       "             │              │  \n",
       "  state_1: ──┼────■─────────┼──\n",
       "             │  ┌─┴─┐┌───┐  │  \n",
       "compare_0: ──┼──┤ X ├┤ X ├──┼──\n",
       "           ┌─┴─┐└─┬─┘└───┘┌─┴─┐\n",
       "     a0_0: ┤ X ├──■───────┤ X ├\n",
       "           └───┘          └───┘"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set x value to estimate the CDF\n",
    "x_eval = 2\n",
    "\n",
    "comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n",
    "comparator.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cdf_circuit(x_eval):\n",
    "    # define the registers for convenience and readability\n",
    "    qr_state = QuantumRegister(u.num_qubits, 'state')\n",
    "    qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')\n",
    "    qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')\n",
    "    qr_obj = QuantumRegister(1, 'objective')\n",
    "    qr_compare = QuantumRegister(1, 'compare')\n",
    "\n",
    "    # define the circuit\n",
    "    state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name='A')\n",
    "\n",
    "    # load the random variable\n",
    "    state_preparation.append(u, qr_state)\n",
    "\n",
    "    # aggregate\n",
    "    state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])\n",
    "\n",
    "    # comparator objective function\n",
    "    comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n",
    "    state_preparation.append(comparator, qr_sum[:] + qr_obj[:] + qr_carry[:])\n",
    "\n",
    "    # uncompute aggregation\n",
    "    state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
    "    \n",
    "    return state_preparation\n",
    "    \n",
    "state_preparation = get_cdf_circuit(x_eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we first use quantum simulation to validate the quantum circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">             ┌───────┐┌────────┐        ┌───────────┐\n",
       "    state_0: ┤0      ├┤0       ├────────┤0          ├\n",
       "             │       ││        │        │           │\n",
       "    state_1: ┤1      ├┤1       ├────────┤1          ├\n",
       "             │  P(X) ││        │        │           │\n",
       "    state_2: ┤2      ├┤2       ├────────┤2          ├\n",
       "             │       ││        │        │           │\n",
       "    state_3: ┤3      ├┤3       ├────────┤3          ├\n",
       "             └───────┘│  adder │┌──────┐│  adder_dg │\n",
       "objective_0: ─────────┤        ├┤2     ├┤           ├\n",
       "                      │        ││      ││           │\n",
       "      sum_0: ─────────┤4       ├┤0     ├┤4          ├\n",
       "                      │        ││  cmp ││           │\n",
       "      sum_1: ─────────┤5       ├┤1     ├┤5          ├\n",
       "                      │        ││      ││           │\n",
       "    carry_0: ─────────┤6       ├┤3     ├┤6          ├\n",
       "                      └────────┘└──────┘└───────────┘</pre>"
      ],
      "text/plain": [
       "             ┌───────┐┌────────┐        ┌───────────┐\n",
       "    state_0: ┤0      ├┤0       ├────────┤0          ├\n",
       "             │       ││        │        │           │\n",
       "    state_1: ┤1      ├┤1       ├────────┤1          ├\n",
       "             │  P(X) ││        │        │           │\n",
       "    state_2: ┤2      ├┤2       ├────────┤2          ├\n",
       "             │       ││        │        │           │\n",
       "    state_3: ┤3      ├┤3       ├────────┤3          ├\n",
       "             └───────┘│  adder │┌──────┐│  adder_dg │\n",
       "objective_0: ─────────┤        ├┤2     ├┤           ├\n",
       "                      │        ││      ││           │\n",
       "      sum_0: ─────────┤4       ├┤0     ├┤4          ├\n",
       "                      │        ││  cmp ││           │\n",
       "      sum_1: ─────────┤5       ├┤1     ├┤5          ├\n",
       "                      │        ││      ││           │\n",
       "    carry_0: ─────────┤6       ├┤3     ├┤6          ├\n",
       "                      └────────┘└──────┘└───────────┘"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_preparation.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Operator CDF(2) = 0.9591\n",
      "Exact    CDF(2) = 0.9591\n"
     ]
    }
   ],
   "source": [
    "# evaluate resulting statevector\n",
    "var_prob = 0\n",
    "for i, a in enumerate(job.result().get_statevector()):\n",
    "    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]\n",
    "    prob = np.abs(a)**2\n",
    "    if prob > 1e-6 and b[0] == '1':\n",
    "        var_prob += prob\n",
    "print('Operator CDF(%s)' % x_eval + ' = %.4f' % var_prob)\n",
    "print('Exact    CDF(%s)' % x_eval + ' = %.4f' % cdf[x_eval])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we run QAE to estimate the CDF for a given $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:    \t0.9591\n",
      "Estimated value:\t0.9595\n",
      "Confidence interval: \t[0.9584, 0.9605]\n"
     ]
    }
   ],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae_cdf = IterativeAmplitudeEstimation(state_preparation=state_preparation,\n",
    "                                      epsilon=epsilon, alpha=alpha,\n",
    "                                      objective_qubits=[len(qr_state)])\n",
    "result_cdf = ae_cdf.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)\n",
    "\n",
    "# print results\n",
    "conf_int = np.array(result_cdf['confidence_interval'])\n",
    "print('Exact value:    \\t%.4f' % cdf[x_eval])\n",
    "print('Estimated value:\\t%.4f' % result_cdf['estimation'])\n",
    "print('Confidence interval: \\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Value at Risk\n",
    "\n",
    "In the following we use a bisection search and QAE to efficiently evaluate the CDF to estimate the value at risk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_ae_for_cdf(x_eval, epsilon=0.01, alpha=0.05, simulator='qasm_simulator'):\n",
    "\n",
    "    # construct amplitude estimation \n",
    "    state_preparation = get_cdf_circuit(x_eval)\n",
    "    ae_var = IterativeAmplitudeEstimation(state_preparation=state_preparation,\n",
    "                                          epsilon=epsilon, alpha=alpha,\n",
    "                                          objective_qubits=[len(qr_state)]) \n",
    "    result_var = ae_var.run(quantum_instance=Aer.get_backend(simulator), shots=100)\n",
    "    \n",
    "    return result_var['estimation']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection_search(objective, target_value, low_level, high_level, low_value=None, high_value=None):\n",
    "    \"\"\"\n",
    "    Determines the smallest level such that the objective value is still larger than the target\n",
    "    :param objective: objective function\n",
    "    :param target: target value\n",
    "    :param low_level: lowest level to be considered\n",
    "    :param high_level: highest level to be considered\n",
    "    :param low_value: value of lowest level (will be evaluated if set to None)\n",
    "    :param high_value: value of highest level (will be evaluated if set to None)\n",
    "    :return: dictionary with level, value, num_eval\n",
    "    \"\"\"\n",
    "\n",
    "    # check whether low and high values are given and evaluated them otherwise\n",
    "    print('--------------------------------------------------------------------')\n",
    "    print('start bisection search for target value %.3f' % target_value)\n",
    "    print('--------------------------------------------------------------------')\n",
    "    num_eval = 0\n",
    "    if low_value is None:\n",
    "        low_value = objective(low_level)\n",
    "        num_eval += 1\n",
    "    if high_value is None:\n",
    "        high_value = objective(high_level)\n",
    "        num_eval += 1    \n",
    "        \n",
    "    # check if low_value already satisfies the condition\n",
    "    if low_value > target_value:\n",
    "        return {'level': low_level, 'value': low_value, 'num_eval': num_eval, 'comment': 'returned low value'}\n",
    "    elif low_value == target_value:\n",
    "        return {'level': low_level, 'value': low_value, 'num_eval': num_eval, 'comment': 'success'}\n",
    "\n",
    "    # check if high_value is above target\n",
    "    if high_value < target_value:\n",
    "        return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'returned low value'}\n",
    "    elif high_value == target_value:\n",
    "        return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'success'}\n",
    "\n",
    "    # perform bisection search until\n",
    "    print('low_level    low_value    level    value    high_level    high_value')\n",
    "    print('--------------------------------------------------------------------')\n",
    "    while high_level - low_level > 1:\n",
    "\n",
    "        level = int(np.round((high_level + low_level) / 2.0))\n",
    "        num_eval += 1\n",
    "        value = objective(level)\n",
    "\n",
    "        print('%2d           %.3f        %2d       %.3f    %2d            %.3f' \\\n",
    "              % (low_level, low_value, level, value, high_level, high_value))\n",
    "\n",
    "        if value >= target_value:\n",
    "            high_level = level\n",
    "            high_value = value\n",
    "        else:\n",
    "            low_level = level\n",
    "            low_value = value\n",
    "\n",
    "    # return high value after bisection search\n",
    "    print('--------------------------------------------------------------------')\n",
    "    print('finished bisection search')\n",
    "    print('--------------------------------------------------------------------')\n",
    "    return {'level': high_level, 'value': high_value, 'num_eval': num_eval, 'comment': 'success'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------------------------\n",
      "start bisection search for target value 0.950\n",
      "--------------------------------------------------------------------\n",
      "low_level    low_value    level    value    high_level    high_value\n",
      "--------------------------------------------------------------------\n",
      "-1           0.000         1       0.751     3            1.000\n",
      " 1           0.751         2       0.960     3            1.000\n",
      "--------------------------------------------------------------------\n",
      "finished bisection search\n",
      "--------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# run bisection search to determine VaR\n",
    "objective = lambda x: run_ae_for_cdf(x)\n",
    "bisection_result = bisection_search(objective, 1-alpha, min(losses)-1, max(losses), low_value=0, high_value=1)\n",
    "var = bisection_result['level']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated Value at Risk:  2\n",
      "Exact Value at Risk:      2\n",
      "Estimated Probability:    0.960\n",
      "Exact Probability:        0.959\n"
     ]
    }
   ],
   "source": [
    "print('Estimated Value at Risk: %2d' % var)\n",
    "print('Exact Value at Risk:     %2d' % exact_var)\n",
    "print('Estimated Probability:    %.3f' % bisection_result['value'])\n",
    "print('Exact Probability:        %.3f' % cdf[exact_var])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional Value at Risk\n",
    "\n",
    "Last, we compute the CVaR, i.e. the expected value of the loss conditional to it being larger than or equal to the VaR.\n",
    "To do so, we evaluate a piecewise linear objective function $f(L)$, dependent on the total loss $L$, that is given by\n",
    "\n",
    "$$ f(L) = \\begin{cases} \n",
    "0 & \\text{if}\\quad L \\leq VaR \\\\\n",
    "L & \\text{if}\\quad L > VaR.\n",
    "\\end{cases} $$\n",
    "\n",
    "To normalize, we have to divide the resulting expected value by the VaR-probability, i.e. $\\mathbb{P}[L \\leq VaR]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">        ┌─────────┐┌──────┐┌─────────┐┌─────────┐\n",
       "q138_0: ┤0        ├┤0     ├┤0        ├┤0        ├\n",
       "        │         ││      ││         ││         │\n",
       "q138_1: ┤1 LinRot ├┤1     ├┤1 LinRot ├┤1        ├\n",
       "        │         ││      ││         ││         │\n",
       "q139_0: ┤2        ├┤  cmp ├┤2        ├┤  cmp_dg ├\n",
       "        └─────────┘│      │└────┬────┘│         │\n",
       "  a4_0: ───────────┤2     ├─────■─────┤2        ├\n",
       "                   │      │           │         │\n",
       "  a4_1: ───────────┤3     ├───────────┤3        ├\n",
       "                   └──────┘           └─────────┘</pre>"
      ],
      "text/plain": [
       "        ┌─────────┐┌──────┐┌─────────┐┌─────────┐\n",
       "q138_0: ┤0        ├┤0     ├┤0        ├┤0        ├\n",
       "        │         ││      ││         ││         │\n",
       "q138_1: ┤1 LinRot ├┤1     ├┤1 LinRot ├┤1        ├\n",
       "        │         ││      ││         ││         │\n",
       "q139_0: ┤2        ├┤  cmp ├┤2        ├┤  cmp_dg ├\n",
       "        └─────────┘│      │└────┬────┘│         │\n",
       "  a4_0: ───────────┤2     ├─────■─────┤2        ├\n",
       "                   │      │           │         │\n",
       "  a4_1: ───────────┤3     ├───────────┤3        ├\n",
       "                   └──────┘           └─────────┘"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define linear objective\n",
    "breakpoints = [0, var]\n",
    "slopes = [0, 1]\n",
    "offsets = [0, 0]  # subtract VaR and add it later to the estimate\n",
    "f_min = 0\n",
    "f_max = 3 - var\n",
    "c_approx = 0.25\n",
    "\n",
    "cvar_objective = LinearAmplitudeFunction(\n",
    "    agg.num_sum_qubits,\n",
    "    slopes, \n",
    "    offsets, \n",
    "    domain=(0, 2**agg.num_sum_qubits - 1),\n",
    "    image=(f_min, f_max),\n",
    "    rescaling_factor=c_approx,\n",
    "    breakpoints=breakpoints\n",
    ")\n",
    "\n",
    "cvar_objective.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7fdc6a68f110>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define the registers for convenience and readability\n",
    "qr_state = QuantumRegister(u.num_qubits, 'state')\n",
    "qr_sum = QuantumRegister(agg.num_sum_qubits, 'sum')\n",
    "qr_carry = QuantumRegister(agg.num_carry_qubits, 'carry')\n",
    "qr_obj = QuantumRegister(1, 'objective')\n",
    "qr_work = QuantumRegister(cvar_objective.num_ancillas - len(qr_carry), 'work')\n",
    "\n",
    "# define the circuit\n",
    "state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, qr_work, name='A')\n",
    "\n",
    "# load the random variable\n",
    "state_preparation.append(u, qr_state)\n",
    "\n",
    "# aggregate\n",
    "state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])\n",
    "\n",
    "# linear objective function\n",
    "state_preparation.append(cvar_objective, qr_sum[:] + qr_obj[:] + qr_carry[:] + qr_work[:])\n",
    "\n",
    "# uncompute aggregation\n",
    "state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we first use quantum simulation to validate the quantum circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = execute(state_preparation, backend=Aer.get_backend('statevector_simulator'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated CVaR: 3.3121\n",
      "Exact CVaR:     3.0000\n"
     ]
    }
   ],
   "source": [
    "# evaluate resulting statevector\n",
    "value = 0\n",
    "for i, a in enumerate(job.result().get_statevector()):\n",
    "    b = ('{0:0%sb}' % (len(qr_state) + 1)).format(i)[-(len(qr_state) + 1):]\n",
    "    am = np.round(np.real(a), decimals=4)\n",
    "    if np.abs(am) > 1e-6 and b[0] == '1':\n",
    "        value += am**2\n",
    "\n",
    "# normalize and add VaR to estimate\n",
    "value = cvar_objective.post_processing(value)\n",
    "d = (1.0 - bisection_result['value'])\n",
    "v = value / d if d != 0 else 0\n",
    "normalized_value = v + var\n",
    "print('Estimated CVaR: %.4f' % normalized_value)\n",
    "print('Exact CVaR:     %.4f' % exact_cvar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we run QAE to estimate the CVaR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae_cvar = IterativeAmplitudeEstimation(state_preparation=state_preparation,\n",
    "                                       epsilon=epsilon, alpha=alpha, \n",
    "                                       objective_qubits=[len(qr_state)],\n",
    "                                       post_processing=cvar_objective.post_processing)\n",
    "result_cvar = ae_cvar.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact CVaR:    \t3.0000\n",
      "Estimated CVaR:\t3.2474\n"
     ]
    }
   ],
   "source": [
    "# print results\n",
    "d = (1.0 - bisection_result['value'])\n",
    "v = result_cvar['estimation'] / d if d != 0 else 0\n",
    "print('Exact CVaR:    \\t%.4f' % exact_cvar)\n",
    "print('Estimated CVaR:\\t%.4f' % (v + var))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-22T01:56:12.651056Z",
     "start_time": "2019-08-22T01:56:12.640412Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>None</td></tr><tr><td>Terra</td><td>0.17.0.dev0+4ada179</td></tr><tr><td>Aer</td><td>0.6.1</td></tr><tr><td>Ignis</td><td>0.5.0.dev0+470d8cc</td></tr><tr><td>Aqua</td><td>0.9.0.dev0+5a88b59</td></tr><tr><td>IBM Q Provider</td><td>0.8.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.7 (default, May  6 2020, 04:59:01) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Tue Oct 20 11:03:45 2020 CEST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
